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ABSTRACT 

We measure the power spectrum, Pp(k,z), of the transmitted flux in the Lya for- 
est using 3035 high redshift quasar spectra from the Sloan Digital Sky Survey. This 
sample is almost two orders of magnitude larger than any previously available data set, 
yielding statistical errors of ~ 0.6% and ~ 0.005 on, respectively, the overall ampli- 
tude and logarithmic slope of Pp(k, z). This unprecedented statistical power requires a 
correspondingly careful analysis of the data and of possible systematic contaminations 
in it. For this purpose we reanalyze the raw spectra to make use of information not 
preserved by the standard pipeline. We investigate the details of the noise in the data, 
resolution of the spectrograph, sky subtraction, quasar continuum, and metal absorp- 
tion. We find that background sources such as metals contribute significantly to the 
total power and have to be subtracted properly. We also find clear evidence for Silll 
correlations with the Lya forest and suggest a simple model to account for this con- 
tribution to the power. While it is likely that our newly developed analysis technique 
does not eliminate all systematic errors in the Pp(k, z) measurement below the level of 
the statistical errors, our tests indicate that any residual systematics in the analysis are 
unlikely to affect the inference of cosmological parameters from z). These results 

should provide an essential ingredient for all future attempts to constrain modeling of 
structure formation, cosmological parameters, and theories for the origin of primordial 
fluctuations. 
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Subject headings: cosmology: data analysis/observations — intergalactic medium — large- 
scale structure of universe — quasars: absorption lines 

1. Introduction 

Although the Lya forest was discovered many decades ago (Lynds 1971), it has only recently 
emerged as one of the prime tracers of the large scale structure in the Universe. The high resolution 
measurements using the Keck HIRES spectrograph (Vogt et al. 1994) have been largely reproduced 
using hydro dynamical simulations (Cen et al. 1994; Zhang et al. 1995; Hernquist et al. 1996; Theuns 
et al. 1998) and semi-analytical models (Gnedin & Hui 1998). The picture that has emerged 
from these studies is one in which the neutral gas responsible for the absorption is in a relatively 
low density, smooth environment, which implies a simple connection between the gas and the 
underlying dark matter. The neutral fraction of the gas is determined by the interplay between 
the recombination rate (which depends on the temperature of the gas) and ionization caused by 
ultraviolet photons. Photoionization heating and expansion cooling cause the gas density and 
temperature to be tightly related, except where mild shocks heat up the gas. This leads to a tight 
relation between the absorption and the gas density. Finally, the gas density is closely related to 
the dark matter density on large scales, while on small scales the effects of thermal broadening 
and Jeans smoothing must be included. In the simplest picture described here, all of the physics 
ingredients are known and can be modeled. The fact that one can trace the fluctuations over a range 
of redshifts (2 < z < 6 using ground based spectrographs) and over a range of scales, which are 
typically smaller than the scales of other tracers, is the main strength of this method. It becomes 
particularly powerful when combined with cosmic microwave background (CMB) anisotropics or 
other tracers that are sensitive to larger scales. Such a combination is sensitive to the shape of the 
primordial spectrum of fluctuations, which is one of the few observationally accessible probes of the 
early universe. These observations are therefore directly testing the models of the early universe 
such as inflation. 

Lya forest observations and constraints on cosmology have been explored by several groups in 
the past. Most of the analyses focused on the power spectrum, Pf(&), of the fluctuations in the 
Lya forest flux, 

5 F (\) = exp[-r(A)]/ (exp(-r)) - 1 , (1) 

where r is the optical depth to Lya absorption. The first such work was by Croft et al. (1998), 
followed by McDonald et al. (2000), Croft et al. (2002), and Kim et al. (2003). These groups were 
limited to a few dozen spectra at most. Recent theoretical analyses, in addition to above, have 
been performed by Gnedin & Hamilton (2002), Zaldarriaga et al. (2001), and Seljak et al. (2003). 
In the latter two of these papers the degeneracy between the amplitude and slope of the primordial 
power spectrum and the normalization of the optical depth-density relation [most sensitive to the 
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intensity of UV background, and typically parametrized in terms of the mean transmitted flux 
fraction, F = (exp(— r))] was emphasized, which leads to a significant expansion of the allowed 
range of cosmological parameters relative to what one would have inferred from the errors on the 
flux power spectrum alone. Seljak et al. (2003) have shown that the current Lya forest constraints 
are consistent with the ACDM model favored by recent CMB data, testing it in a regime of redshift 
and length scale not probed by other measurements, but that within the ACDM framework they 
do not add much leverage on parameter values beyond that afforded by the CMB data alone. 

An important practical implication of the theoretical breakthroughs of the 1990s is that large 
scale structure in the Lya forest can be effectively studied with moderate resolution spectra. Once 
the spectrum is modeled as a continuous phenomenon rather than a collection of discrete lines, there 
is no need to resolve every feature. Some of the studies cited above use high resolution (~ 0.08A) 
spectra, some use moderate resolution (~ 1 — 3A) spectra, and some use a mix of the two. 

The goal of this paper is to present a new measurement of the Lya forest power spectrum, based 
on ~ 3000 Sloan Digital Sky Survey (SDSS; York et al. (2000)) spectra that probe the Lya forest at 
a resolution R ~ 2000 (~ 2. 5 A FWHM). This sample is almost two orders of magnitude larger than 
anything that was available before. As such it greatly increases the statistical power of the Lya 
forest, making it comparable to the CMB from WMAP. At the same time, the required tolerance 
of systematic errors also increases by the same amount. This requires a careful investigation of 
all of the sources of systematic errors, and a large portion of this paper is devoted to the issue of 
possible systematics in the data and their influence on the parameters of interest. We also discuss 
how the analysis we perform and results we obtain differ from what can be done using the standard 
spectral pipeline outputs in the public SDSS data. In part because of the practicalities of work in 
a large, multi-institutional collaboration, and in part because of the importance of obtaining an 
accurate measurement with well understood statistical and systematic errors, the Lya forest power 
spectrum has been pursued by two independent groups within the SDSS, one led by P. McDonald 
and U. Seljak, and the other by L. Hui and A. Lidz. The methods employed are different and have 
been developed independently. Results of the alternative analysis will be presented elsewhere (Hui 
et al., in preparation). 

We only present the observational measurement of the SDSS Lya forest power spectrum in the 
current paper. Independent of any theoretical interpretation, this basic result should be robust on 
the scales for which we give results, 0.0013 ( kms -1 ) -1 < k < 0.02 ( kms -1 ) -1 , where k = 2ir/\ if A 
is the wavelength of a Fourier mode (not to be confused with spectral wavelength) , here measured in 
kms -1 (note that throughout the paper we frequently use velocity in place of observed wavelength, 
with the understanding that all that enters into our calculations are velocity differences between 
pixels of measured spectra, defined by Av = c A ln(A) - we do not measure power on scales large 
enough for the imperfections in this expression to become relevant). The choice of /c-range is 
determined by the continuum fluctuations on the low end and spectral resolution at the high end. 
We note also that the useful range is limited not only by these uncertainties, which are related 
to the data analysis, but also by the uncertainties in the theoretical modeling and/or additional 



-4 - 



astrophysical effects. We will address these latter issues in more detail in a separate publication. 
However, we do not completely decouple the theory from the data analysis. For example, when 
discussing the importance of systematic errors it is useful to understand how they would affect 
cosmological results like the slope and amplitude of the matter power spectrum, so much of our 
discussion of systematics is devoted to this issue. 

The common usage of the term Lya forest is to describe the Lya absorption by neutral hydrogen 
in the relatively low density bulk of the IGM. In this paper we include damped-Lya systems (DLAs) 
in the definition of the "forest" , so it includes all HI-Lya absorption. We could try to remove DLAs 
before measuring Pp(k, z), because they are more difficult to simulate than the lower optical depth 
absorption; however, we believe the advantage of removal is illusory. If the DLAs were located 
randomly within the IGM (which they certainly are not completely), it would be simple to include 
them in the theory using their known column density distribution. If they are not located randomly, 
the regions obscured by DLAs in the spectra are special, so the effect of removing the DLAs still 
must be understood using simulations. We leave the handling of the effects of DLAs as a problem 
for the theory, which we will address elsewhere. 

Absorption by metals is also difficult to simulate accurately, so we would like to remove this 
contribution to Pp{k,z). This is relatively easy to do for transitions with wavelength A > 1300 
A, but it is basically impossible for transitions with A < X a , because the metal features always 
appear mixed with HI-Lya. We will subtract the power measured in the rest wavelength range 
1268 A < A rcs t < 1380 A from our measurement of the power in the forest, which removes the 
effect of transitions with longer wavelength, but we leave shorter wavelength transitions as part of 
the forest. The only significant contaminant of this kind that we can identify is Silll absorption 
at 1206. 50A, and we develop a simple and effective way to account for this in the theory. We refer 
to our final background-subtracted power spectrum as Pp(k,z), and use P\ 1 \ 2 (k, z) for the raw 
power measured in the interval Ai < A rcs t < A2. We are using the range 1041 A < A rcst < 1185 A 
for the Lya forest. 

The outline of this paper is as follows. In §2, we describe the selection of our data set and 
the preparation of spectra for the measurement of Pp(k, z). In §3 we describe the method used to 
measure the power spectrum and estimate the error bars, test the procedure, and give the basic 
results. We perform consistency checks on the results and discuss systematic errors in §4, which is 
followed by a brief recipe for using our results in §5, and conclusions in §6. 

2. Data Selection and Preparation 

We describe the sample of quasar spectra that we use in §2.1. In §2.2 we explain how we 
remove broad absorption line (BAL) quasars from the sample. In §2.3 we explain how we combine 
spectra from different exposures for the same quasar and use the differences between exposures to 
understand the noise in the data. We discuss the resolution of the spectra in §2.4. Finally, in §2.5 
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we describe how we divide each spectrum by an estimate of the quasar continuum, the expected 
mean absorption level in the spectrum, and a spectral calibration vector (see below), to produce 
the vectors of transmission-fluctuation estimates, Sf, for each quasar, from which we will measure 
P F (k,z). 

2.1. SDSS Observations and Sample Selection 

The Sloan Digital Sky Survey (York et al. 2000) uses a drift-scanning imaging camera (Gunn 
et al. 1998) and a 640 fiber double spectrograph on a dedicated 2.5 m telescope. It is an ongoing 
survey to image 10,000 sq. deg. of the sky in the SDSS ugriz AB magnitude system (Fukugita 
et al. 1996; Stoughton et al. 2002) and to obtain spectra for ~ 10 6 galaxies and ~ 10 5 quasars. The 
astrometric calibration is good to better than O'/l rms per coordinate (Pier et al. 2003), and the 
photometric calibration is accurate to 3% or better (Hogg et al. 2001; Smith et al. 2002). The data 
sample used in this paper was compiled in Summer 2002 and is a combination of data releases one 
(Abazajian et al. 2003) and two (Abazajian et al. 2004). 

About 13% of the spectroscopic survey targets are quasar candidates selected based on their 
colors (Richards et al. 2002). The magnitude limit for UV-excess objects is i = 19.1, while additional 
high-redshift candidates (z > 3) are targeted to i = 20.2. Fibers are allocated according to a tiling 
algorithm (Blanton et al. 2003), with the galaxy sample and the quasar sample being the top 
priorities. The remaining 8% of fibers serve for calibration purposes. 

SDSS spectra are obtained using plates holding 640 fibers, each of which subtends 3" on the 
sky; the spectra cover 3800 - 9200A. The pixel width is a slowly varying function of wavelength, 
but is typically ~ 70kms _1 . The resolution also varies, but is typically also ~ 70kms _1 rms (i.e., 
the resolution is 1800 < R < 2100 and there are ~ 2.4 pixels per FWHM resolution element). 
All quasars have multiple spectra, usually taken one after the other (timescales of a fraction of an 
hour), so the quasar variability can be ignored (in the opposite case it would act as an additional 
source of noise). The co-added spectra in the official SDSS release use local spline interpolation 
onto a uniform grid of pixels of width Alog 10 (A) = 0.0001, and do not guarantee the noise to be 
uncorrelated. We therefore redo this step starting from the individual exposures. This is discussed 
in more detail below. Spectral flux errors per pixel in most cases are about 1 x 10~ 17 erg s _1 cm~ 2 
A -1 . Redshifts are automatically assigned by the SDSS spectral classification algorithm, which is 
based on \ 2 fitting of templates to each spectrum (Schlegel et al., in preparation). 

We limit ourselves to quasars with redshift z q > 2.3 when measuring the power in the Lya 
forest region of spectra, so that each spectrum contains a significant stretch of the Lya forest above 
the detector cutoff at 3800 A (which corresponds to Lya absorption at z = 2.12). We use the sample 
compiled in Summer 2002, cut down to 3035 spectra by eliminating some plates of questionable 
quality, some spectra where two different redshift estimation codes disagree, and some BAL quasars 
(see below). Figure 1 shows the redshift distribution of the data. The dashed, red histogram shows 



- 6 - 



1.5x10 



10 5 - 



m 

O 
X 

• r 1 



5xl0 4 



- 300 




500 



400 



C 
GO 
"I 



200 



100 



Fig. 1. — The distribution of the spectral pixels used to probe the Lya forest (black, solid histogram; 
scale on left axis), and the redshift distribution of our primary sample of 3035 quasars (red, dotted 
histogram; right axis). Note the gap at z ~ 2.7 in the quasar redshift distribution, caused by a 
class of stars being indistinguishable from quasars in the SDSS photometry (Richards et al. 2002). 
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the distribution of quasar redshifts. The solid, black histogram shows the distribution of pixels in 
the range 1041 < A rest < H85A. Note that there is a gap in the quasar redshifts around z ~ 2.7, 
which is due to the stellar locus crossing the quasar locus in the 5-color SDSS photometry (Richards 
et al. 2002). Figure 2 shows an example SDSS spectrum of a z = 3.7 quasar. This spectrum is 
unusual in that most have lower S/N, and most quasars are at lower redshift. 

We employ an additional sample of ~ 8000 spectra with z q < 2.3, so that we can study the 
full observed wavelength range, 3800 A < A < 9200 A, outside the confusion of the Lya forest. As 
we discuss in §3.4, we compute a non-negligible background power term (probably mostly metal 
absorption), by measuring the power in the wavelength range 1268 A < A rcst < 1380 A. Using only 
the primary sample, we would not be able to compute this term for observed wavelengths below 
~ 4400 A. 

We remove several wavelength regions from our analysis because of calibration problems: A < 
3800 A, 5575 A < A < 5583 A, 5888 A < A < 5894 A, 6296 A < A < 6308 A, and 6862 A < A < 6871 
A (the last two have no direct effect on the results we present). Most of these problems are due to 
strong sky lines. 

2.2. BAL Removal 

Our sample was initially examined by eye, and the most extreme broad absorption line (BAL) 
quasars were removed (see Hall et al. (2002) for a discussion of BALs). When we first measured 
the background power in the region 1409 < A rcst < 1523A, we found that the most extreme 
outliers in power were still obvious BALs (this was not true of the Lya forest region). To test 
the importance of these systems to our Lya forest power measurement, we removed a further 147 
quasars using the following automated method: Each spectrum is smoothed by a Gaussian with 
rms width 280 km s -1 . The continuum within the region 1420 < A rcst < 1535A is redefined by 
dividing by the mean flux-to-continuum ratio in the region. A quasar is identified as a BAL quasar 
if the region 1420 < A rcst < 1535A contains a 2000kms _1 long continuous set of pixels that all 
fall more than 20% above or below our estimated continuum (we initially identified wide regions 
with flux above the continuum out of simple curiosity, but found that these are in practice almost 
always obvious BAL quasars where the continuum has been biased low by the BAL feature). We 
iterate the continuum redefinition twice, computing the new mean after throwing out pixels more 
than 20% below the previous mean, but this makes almost no difference to the results. Note that 
the 280 km s -1 smoothing was applied to allow easier identification of BALs in noisy spectra. As 
we show below, this BAL cut makes essentially no difference to our Pp(k,z) result, although it 
does have a noticeable effect on the power measurement in the region 1409 < A rcs t < 1523A. 
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Fig. 2. — Example spectrum of a z = 3.7 quasar with unusually high S/N. The regions we use to 
measure the Lya forest power and background power are indicated by vertical dotted lines, along 
with a couple of alternate regions that we will discuss (note that the background and Lya forest 
observed in the same quasar spectrum correspond to different redshifts). 
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2.3. Combining Exposures and Calibrating the Noise 

SDSS obtains multiple (at least 3) exposures for each quasar. We combine the individual 
exposure spectra to produce a single spectrum, using a nearest-grid-point method that produces 
uncorrelated noise and a reasonably well-defined sampling window. For each pixel we record es- 
timates of wavelength, quasar flux, resolution, sky flux, read-noise, and two different total noise 
estimates. The first noise estimate, which we will call simply a v (p for pipeline), is computed 
using the error array given for the exposure spectra by the spectral reduction pipeline. The second 
noise estimate, which we call o c (c for component), is computed by summing the read- noise and 
the noise implied by estimates of the number of photons corresponding to the quasar flux and sky 
flux. The two noise estimates generally do not agree, but this is not a problem for us because 
we ultimately recalibrate the noise (next). Finally, we record \ 2 l v f° r each pixel, computed by 
treating the determination of the combined flux value for each pixel as a one parameter fit to the 
measurements given by the different exposures. Examples of the more important of these quantities 
in Lya forest regions are shown in Figure 3. For comparison to the sky and quasar flux levels, we 
have converted the Gaussian read-noise into the flux of photons that would contribute the same 
noise variance. Several elements of Figure 3 (e.g., the estimation of the quasar continuum and a w ) 
will be described later in this paper. 

The noise estimate from the standard SDSS pipeline is only approximate. The accuracy of 
the noise estimate required for our purpose is much higher than anticipated when the pipeline was 
developed. For this reason we use the differences between single-exposure spectra for the same 
quasar to determine the noise properties of the data. We construct difference spectra by combining 
the flux-calibrated exposures with alternating sign for each exposure, i.e., we use exactly the same 
procedure that we normally use to produce combined spectra from the exposures, except half of the 
exposures are subtracted instead of added, so the mean result is zero (we drop the last exposure 
when there are an odd number - this is not the most efficient method possible, but we do not 
need it to be). The result is a direct measure of the exposure-to-exposure changes. We measure 
the power spectrum of these difference spectra using the method described in §3.1, including noise 
subtraction based on the pipeline noise estimates for the pixels. The result is shown in Figure 4 
(points with error bars). We obtain a clear detection of power, where there should be none if the 
spectra differ only by the noise estimate from the pipeline which is being subtracted. If we assume 
that the noise has been underestimated by a constant factor, and fit for that factor using the error 
covariance matrix estimated by bootstrap resampling, we find a decent fit: \ 2 = 141.6 for 107 
degrees of freedom (formally, this fit is not good because x 2 1S unlikely to be this high by chance). 
This fit uses our usual points in 0.0013 (kms -1 )" 1 < k < 0.02 (kms -1 )" 1 . The best fit value of 
the excess noise contribution is 16.1 ± 0.4% of the original noise estimate, indicating that the rms 
noise was underestimated by 8%. The best fit and goodness of fit do not change if we add points on 
larger scales. The quality of the fit begins to degrade as we add points with larger k, but the best fit 
value changes by only 1% (in power) out to the Nyquist frequency of the data. Of course, we have 
no reason to expect a single redshift-independent factor to describe the relation between the true 
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Fig. 3. — Examples of the chunks of spectra used to measure power, with (a,b) showing quasars 
at z q =(3.24,2.45) over the rest wavelength range 1113 A < A rcst < 1185 A, and (c) showing a 
quasar at z q = 3.30 over the rest wavelength range 1041 A < A rest < 1113 A. Top panel: quasar 
flux (solid black line), sky flux (dotted blue line), our continuum estimate (red short-dashed line), 
and the read-noise as an equivalent photon flux (green long-dashed line). Middle panel: S/N level 
shown as a ratio of our continuum to the different rms noise levels (see text), a w (black solid line), 
a p , (blue dotted line), and a c (red dashed line). Bottom panel: Calibration correction vector, S 
(blue dotted line), rms resolution in units of 100 km/s (red dashed line), and evolution of the mean 
transmission fraction, F(z) (black solid line). The perfect degeneracy in our analysis between the 
overall normalization of the continuum and F(z) has been broken arbitrarily, so only the evolution 
of F(z) is meaningful (see text). 
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Fig. 4. — The points show the measured power in difference spectra, created by subtracting separate 
exposures for the same quasar. Noise power has been subtracted based on the standard pipeline 
noise estimates for each exposure. The lines show 16% of the subtracted noise term. The different 
colors, lines, and symbols identify redshift bins, in a pattern that we will use repeatedly throughout 
the paper. From bottom to top — z=2.2: black, solid line, open square; z=2.4: blue, dotted line, 
4-point star (cross); z=2.6: cyan, dashed line, filled square; z=2.8: green, long-dashed line, open 
triangle; z=3.0: magenta, dot-dashed line, 3-point star; z=3.2: red, dot-long-dashed line, filled 
triangle; z=3.4: black, thin solid line, open pentagon; z=3.6: blue, thin dotted line, 5-point star; 
z=3.8: cyan, thin dashed line, filled pentagon. The different redshifts have been shifted vertically 
by arbitrary amounts on this logarithmic plot. 
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and pipeline noise, so the formally bad x 2 1S n °t a fundamental problem. We check for systematic 
change with redshift by allowing a power law dependence, -Presidual noise ^ 

[3.75/(1 + z)] d , but find 

no significant detection (d = 0.07 ± 0.20). Our final method will effectively account for evolution 
anyway, as described below. 

A k dependence different than expected for white noise could be a problem for us, so we check 
for this by fitting for a power law dependence, -Presidual noise °c (k/k p ) b [with k p = 0.0074 (kms -1 ) -1 ], 
finding b = —0.111 ± 0.025, a significant detection (x 2 is now a reasonable 123.3 for 106 dof). Al- 
lowing a running of the power law, Presidual noise oc {k/k p ) h+l / 2 c M fc /M, does not improve the fit 
(c = —0.046 ± 0.066). The slope we find corresponds to a ~ 20% change in 16% of the noise power 
at the extremes of our k range, i.e., only ~ 3% of the total noise power, which is a relatively small 
fraction of the Lya forest power except at the highest k (see Figure 1 1 below) . We henceforth as- 
sume that the extra noise is proportional to k~°' in rather than white (this makes < 1% difference 
in the final results except for the one highest k, lowest z point where the difference is 2%). 

How accurate is this noise estimate based on differences between exposures? Our difference 
spectra will contain a component of the Lya forest power if the calibration between exposures 
is not perfect. The power in this term would be suppressed relative to the Lya forest power by 
the fractional calibration error squared, so it would be very small unless the exposure-to-exposure 
calibration errors were quite large. The fact that a simple one parameter extra-noise model fits 
reasonably well, in the face of variation in redshift, noise amplitude, and k, argues against cali- 
bration errors being a big problem. More convincingly, we measure nearly the same excess noise 
contribution (14.2 ± 0.5%) and slope (b = -0.135 ± 0.028) in the region 1268 A < A rC st < 1380 A 
as we do in the Lya forest. This argues against any connection to leaking Lya forest power. Note 
that the effective absolute level of noise in the 1268 A < A rC st < 1380 A region is about half that in 
the Lya forest region, so this test shows that the fraction of extra noise does not depend strongly 
on the noise level itself. 

Pixels in different exposures are not perfectly aligned, and misalignment can allow Lya forest 
power to leak into our difference spectra. To test this alternative explanation for the apparent 
excess noise in the spectra, we split the spectra into two groups with approximately equal weight, 
based on the rms misalignment in the forest region (the alignment is known from the wavelength 
calibration of the exposures, which is thought to be practically perfect). We find the same excess 
noise power in both the poorly aligned group (16.1 ± 0.6%, b = —0.086 ± 0.036) and the better 
aligned group (15.3 ± 0.6%, b = —0.123 ± 0.036), suggesting that the excess power is not due to 
misalignment. Furthermore, the presence of a similar level of excess noise power outside of the forest 
region again argues against leakage. We therefore believe that our noise estimate is considerably 
more accurate than the noise estimate from the SDSS pipeline. 

In our initial power spectrum analysis we multiplied the noise-power estimated from the 
pipeline errors by the factor 1.16 for all spectra; however, when we split the data based on the 
mean value of x 2 l v f° r the exposure combination (see §4.4) we found that the large and small 
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X 2 /v subsamples disagreed significantly on the Pp(fc,z) results. We eliminated this problem by- 
estimating the noise-correction factor individually for each spectrum, by fitting to the power in 
the difference spectrum for that quasar. The mean extra power from these fits is still close to 
16%, but there is considerable scatter. When we use these individual estimates, the correlation 
between measured Pp(k,z) and x 2 / v disappears, i.e, the mean value of x 2 /^ f° r a spectrum's 
exposure combination was a good indicator of the amount by which the noise in each spectrum 
was misestimated. Note that there are statistical errors in these noise estimates for each spectrum, 
of the same order as the error for which we are trying to correct; however, there is no systematic 
bias associated with these errors, and the random error they contribute is automatically included 
in our final bootstrap errors. In fact, including the spectrum-by-spectrum noise estimate reduces 
the bootstrap errors slightly on small scales, verifying that these estimates are on average more 
accurate than the original noise estimates. It is not known why the noise is misestimated by the 
standard pipeline. Tests at this level have not been done before. 

Our final data product will be a measurement of Pp(k, z)binned in k and z, i.e., a matrix PF,ij 
where i labels bins with Zi and j labels bins with kj. We will also give the noise power that was 
subtracted, PN,ij, hi the same bins. We suggest allowing a 5% rms freedom in the noise amplitude 
in each z bin when performing model fits, i.e., for each bin subtract fiPN,ij from PF,ij, and add 
(/i/0.05) 2 to x 2 ■ This is probably overly conservative for any one bin, but implies a combined 
freedom ±0.05/3 (for 9 bins) on an overall noise misestimation. This seems prudent, even though 
it is not really required by any test we have performed. 

2.4. Accuracy of the Resolution 

The resolution of the SDSS spectra is estimated using lines from calibration lamps mounted 
on the telescope structure. Shifts of the detector pixel grid relative to a fixed observed wavelength 
frame during an exposure, which we will call flexure, are expected to be the dominant source of 
error in this spectral resolution estimate. We tried estimating the rate of shifting for each pixel by 
differencing the wavelength calibrations of adjacent exposures (this calibration is determined very 
precisely for each exposure using the positions of sky lines). The implied extra smoothing only 
changes the power by ~ 2% at our highest k bin. 

The strong sky line at 5577 A provides a good opportunity to measure the resolution more 
directly (note that the spectral wavelengths are in vacuum, and heliocentric, so this and other sky 
lines generally appear shifted from their standard wavelength). We measure the power spectrum in 
~ 3000 sky spectra in the range 5560A < A < 5598A. If the sky line has negligible width and the 
smoothing has a Gaussian shape with rms width R, the power spectrum should be proportional to 
W 2 (k,R,l) = exp[-(/c-R) 2 ][sin(/c 1/2) /{k 1/2)] 2 , where I is the pixel width (the pixelization effect 
is sub dominant but not negligible). In Figure 5 we show the measured power averaged over all 
the sky spectra after dividing each individual measurement by W 2 (k, R, I), where R and I are the 
local values (they are to a good approximation constant over the range we are looking at), and also 
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dividing each measurement by the value at a low k where the resolution should not have any effect. 
The result is remarkably close to unity, indicating that the estimated resolution is an accurate 
representation of the true resolution. What are the small wiggles? Figure 6 shows an example of 
the region we Fourier transform to measure the power. We believe the small features to the sides of 
the main line are OH lines at 5564 and 5589 A (Slanger et al. 2003). We test this explanation for 
the wiggles by constructing mock sky spectra that simply have a delta function at 5579A and two 
more with 0.003 times the main line's amplitude at 5566 and 5591 A, convolved with the resolution 
and pixelization (0.003 was chosen to give the best fit to the wiggles). The red, dotted line in 
Figure 5 shows the same resolution test using the mocks. We see that the wiggles are essentially 
perfectly reproduced. In conclusion: the resolution profile appears to be perfectly Gaussian, with 
exactly the width expected from the given resolution. There is apparently no room for even a 2% 
level effect from flexure. We are prevented from performing the same kind of measurement using 
other sky lines by similar features which are always much larger relative to the central line. 

We suggest that fits to Pp(/c,z) include a multiplicative uncertainty on the overall power, of 
the form exp(afc 2 ), where a is a single parameter in the fit subject to the rms constraint a a = 
(7kms -1 ) 2 . This allows for a ~ 2% change in the smoothing kernel at our highest k, similar to our 
estimate of the error from flexure. This error estimate is somewhat arbitrary, but the evidence we 
have presented suggests that it should be smaller, so our estimate is conservative. 

Note that this resolution test, and the noise calibration, cannot be used directly with the 
standard pipeline spectra, where the exposures are combined in a different way. The reader may 
be confused at this point about how our spectra differ from the standard publicly available set, so 
we give the following summary: 

• Our nearest-grid-point combination of the exposures produces uncorrelated noise in pixels (to 
the extent that the noise in the exposures was uncorrelated, which is expected from the way 
they are extracted), while the standard pipeline uses a local splining procedure which does a 
good but not perfect job of preventing noise correlation. 

• When combining exposures we record the effect of different pixel sizes, misalignment of the 
pixels, and flexure of the detector during exposures, which can influence the effective resolu- 
tion. 

• We record the contribution of quasar flux, sky flux, and read-noise to the total noise in 
each pixel. Knowing the contribution from quasar flux is important if pixel-by-pixel noise 
weighting is to be used, because the correlation between flux level and noise amplitude can 
lead to biases (see the end of §2.5). 

• We correct for the bias in the exposure combination associated with cross-correlation be- 
tween the noise variance level in exposure pixels and the quasar flux in the pixel (a different 
incarnation of the problem alluded to in the previous point). 
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Fig. 5. — Resolution test. The solid, black line with error bars shows the power measured in ~ 3000 
sky spectra in the range 5560A < A < 5598A (dominated by the strong sky line at 5577A) divided 
by the asymptotic small k power and by the estimated resolution/pixelization kernel W 2 (k) for 
each spectrum. If the resolution estimate was perfect, and the sky line was narrow and the only 
flux present, this division would give exactly 1. The large error bars are the spectrum-to-spectrum 
variation, the small ones are the error on the mean. The blue, long-dashed line shows the power 
not divided by W 2 (k), i.e., basically an averaged version of W 2 (k), which drops to ~ 0.25 by 
k = 0.02( kms" 1 )" 1 . The red, dotted line shows the result of our test for mock spectra constructed 
with a Gaussian at 5579A and two more at 5566 and 5591 A with 0.003 times its amplitude, 
representing OH lines. The green, short-dashed line shows exp[(fc 7kms -1 ) 2 ]. The vertical, cyan, 
dotted lines bound the k region in which we will present Lya forest results, while the horizontal, 
cyan, dotted line just guides the eye to 1. 



-18- 



0.03 



w 

cd 



in 



0.02 - 



0.01 




Q I I I I i I i I I I I I I I I I I I L 

5560 5570 5580 5590 

A[A] 



Fig. 6. — Example of the sky flux near the sky line at vacuum wavelength 5579A, relative to the 
peak of the line (one spectrum only). 
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• The noise is recalibrated for each spectrum by differencing the exposures. The noise variance 
in the standard pipeline exposure spectra is underestimated by on average 16%, on top of any 
error related to the exposure combination, and the power measured in the difference spectra 
is slightly tilted relative to white noise. 

The last point is the most important. 

2.5. Determination of the Continuum and Mean Absorption Level 

Our goal is to measure the power spectrum of the fluctuations in the transmitted flux fraction 
through the IGM, <5p(A) = F(X)/F — 1, where F(X) = exp[— r(A)] and r(A) is the Lya forest optical 
depth (as defined in §1). However, the spectrum of each quasar is the product of F and the quasar 
continuum (note that we use "continuum" to refer to all the flux emitted by the quasar, including 
emission lines), further complicated by errors in the detector calibration and absorption by longer 
wavelength transitions. The details of the procedure we use to separate these contributions will be 
presented elsewhere, here we give the basic idea and key results that are relevant for the flux power 
spectrum determination. 

We use an iterative procedure to determine the components of the data model 

f = A q C(Xl cst ) (1 + Sh) SiX') (1 + Si) F(z*) (1 + 4) + n* , (2) 

where f % is the raw flux in pixel i, n % is the noise, A q is the overall normalization of the quasar 
spectrum, ^(A^t) is the mean quasar continuum shape, 5c are fluctuations around the mean 
continuum, S(X) is a mean generalized calibration vector (this includes wavelength dependent 
calibration errors in the SDSS spectra, but also the mean absorption by metal lines with resonance 
wavelength A > 1300A), 5,s are fluctuations around S, such as individual metal lines or variable 
calibration errors, F{z) is the mean Lya forest absorption at a given redshift, and 5p are the 
fluctuations in Lya forest absorption. Note that here, as most places in this paper, z l = X 1 / X a — 1 
is the redshift of gas that would produce Lya absorption in the pixel, not the redshift of the quasar. 
Briefly, we determine A q for each spectrum, the global functions C{X Tes t), S(X), and F(z), and a 
set of principal component analysis (PCA) eigenvectors that describe 5c by, in turn, computing 
each component of the model from all the spectra while holding all the others fixed. E.g., F{z) is 
estimated in bins of z by averaging f' l /A q C(Xl est ) S(X l ) (1 + S' c ) over all the pixels in the Lya 
forest that fall in each bin. We separate S(X) from F(z) by measuring S(X) in the rest wavelength 
range 1268 A < A rcst < 1380 A, i.e., outside the Lya forest. A few iterations suffice to determine 
all the components of the model independently. Three examples of the results are shown in Figure 
3. The full details of this procedure will be presented in a separate paper focused on a precise 
determination of F(z). 

In preparation for measuring the power spectrum, we divide each quasar spectrum by our 
estimate of A q C(X rcst ) (1 + 5c) S(X) F(z). The power in the 5(A) and F(z) terms is completely 
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negligible (always < 0.5% of Pf(&, z) and usually much less). 5c are represented for each quasar by 
N PCA eigenvectors. We have tried several different values for N ranging from to 13, and find that 
the power spectrum results depend slightly (but not critically) on the value, as we discuss below. 
For our final results we use N = 0, i.e., we only divide by a mean continuum, although we will also 
show how using N = 13 affects the cosmological fit results. We do not use the PCA continua because 
we are not satisfied with their robustness, and division by them frequently actually increases the 
resulting power slightly. This may indicate that the error introduced by allowing additional freedom 
in the continuum is larger than the continuum fluctuations that we are trying to remove. Our study 
of PCA continua was primarily aimed at determining F{z) rather than the power spectrum, so we 
cannot be certain that the PCA method could not be used productively in a power spectrum 
measurement if it was more carefully optimized for that application. Because we know that our 
continuum estimate (which involves an extrapolation from outside the Lya forest region) is not 
perfect within the Lya forest, we further divide each chunk of spectrum that will be used to make 
a power spectrum estimate by its mean (optimally computed considering both observational noise 
and absorption variance). We call our resulting observed data vector 5f = 5f + 5 s + 5 n , where 5 n is 
the normalized noise fluctuation and we are ignoring the cross-terms between 5f, 5c, and 5s- As we 
describe in detail below, 5s is treated as a random noise background and its statistical properties are 
determined by measuring the power spectrum in the rest wavelength range 1268 < A rcst < 1380A, 
where 5f = (and F = 1). 

A small but non-negligible detail of our procedure is hidden within our description of the 
normalization of the spectra. When we estimate the mean to divide by, we weight the computation 
optimally using the covariance matrix, Cij, of the pixels (C is discussed in more detail in our 
explanation of the power spectrum measurement below). C includes Lya forest fluctuations and 
measurement noise. We do not use our best estimate of the measurement noise directly for the 
weighting, because the noise variance estimate is correlated with the measured flux in the pixel, 
which leads to a bias: the mean is underestimated because lower flux pixels have lower noise. The 
original noise estimate is o 2 p = 7 (f quaS ar + fsky) + a readnoise^ where f qua sar is the flux from the 
quasar, f s k y is the flux from the sky, and 7 accounts for the conversion between the units of flux 
and photons (this description is slightly idealized since the reduction of two-dimensional CCD data 
to a spectrum introduces some complications). To remove the correlation between flux and noise, 
we subtract -ff quaS ar horn o 2 v and add 7 (f quaS ar}, where (f qU asar) = \ C{\ x ^) (l + 6 c ) S(X) F(z). 
We call the final result a w (w for weight; see Figure 3 for some comparisons of the noise estimates). 
The estimate of 7 we have from the spectral reduction pipeline is not perfect, so our replacement of 
the correlated part of the noise amplitude is imperfect. We make a final, very small, correction to 
the mean estimation based on a direct computation of the cross-correlation between the flux and 
noise amplitude. We use the same decorrelated noise amplitudes for weighting the power spectrum 
extraction (discussed below); however, the bias is completely insignificant in that case (i.e., Pp(k, z) 
computed using the original noise estimates for weighting is practically identical). 
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3. Power Spectrum Determination 



The high precision of the Pp(k,z) measurement obtainable using the SDSS data sample re- 
quires unprecedented (in this field) care in the design and testing of the procedure used to produce 
it. We describe the basic method that we use to extract a power spectrum and estimate the errors 
in §3.1. In §3.2 we present the test of the full method as implemented in our code, using mock data 
sets. In §3.3 we give the raw result for the measurement of power in the Lya forest region. 

We aim to measure the power spectrum of Sf, representing the correlation of fluctuations in 
the Lya forest absorption only; however, the covariance matrix of the data vector Sf is 



(the three components of 5f as we have defined it should be uncorrelated). The noise term in 
equation (3) is relatively easy to compute and subtract. We estimate and subtract most of (5s5s T ) 
by defining 



where z is always defined by z = X/X a — 1, so that we are subtracting power measured in the 
same observed wavelength ranges, not the same quasar spectrum (we remind the reader that we 
have defined P\ lt \ 2 to mean the power measured in the range Ai < A rcst < A2). As it appears 
in Pi268,i38o(k, z), z is the redshift of gas that would produce Lya absorption in this part of 
the quasar spectrum, if it was not at a higher redshift than the quasar, i.e., z is really just an 
indicator of observed wavelength. The subtraction in equation (4) will completely remove the 
power due to transitions with A > 1380 A, including SilV (a doublet absorbing at rest wavelengths 
1393.75 and 1402.77 A) and CIV (another doublet at 1548.20 and 1550.78 A). Note that this 
subtraction of metal power is exact, not an approximation [except for the approximation that 
(1 + #f)(1 + 5s) — 1 + Sf + Ss], because we are determining the metal power in exactly the same 
observed wavelength range as the Lya forest power from which it is being subtracted, i.e., the 
same gas, at the same redshift, is doing the absorbing both inside and outside the forest, so the 
absorption will have identical statistical properties. This background subtraction will also remove 
any strictly observed-wavelength-dependent power introduced by the detector, such as spectrum to 
spectrum variations in the calibration of the detector. We implement it in §3.4. 




(3) 



PF(k,z) = Pl041,1185(fc,2) - -Pl268,1380(M) , 



(4) 



3.1. 



Core Method 



In this subsection we describe our method for extracting the power spectrum, Px(k,z), from 
any selected rest wavelength range X (§3.1.1), and estimating its statistical uncertainty (§3.1.2). 
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3.1.1. Band-Power Estimation 

We estimate Px(k) using the quadratic estimation method, which is essentially a fast iterative 
implementation of the maximum likelihood estimator (we follow the expressions as given in Seljak 
(1998)). This method is optimal for a Gaussian probability distribution. While the power spectrum 
estimates are not Gaussian distributed, the deviations are small, as shown below. We measure the 
power in flat bands with edges given by log 10 (fcj) = —4.2 + 0.1 i where i ranges from to 30 (to 
produce 29 bands), although we will not give results for some of the large and small-scale bands 
when we think they are unreliable. Defining Sf = 5x + b~n, where 5x are the fluctuations we are 
measuring (e.g., 5x = 8f + <5s within the forest) and 5 n are the normalized noise fluctuations, a 
band-power estimate, P k , for each chunk of spectrum is given by 

p k = \Y. F kk^f T C^Qk'C- l 5f - M , (5) 

k' 

where C = (5 f 5 f T ) = S + N,S= (5 X 5 X T ), N = (5 n 5 n T ), Q k = 8S/dP k , 

F kk < = ^{C^QkC^Qv) (6) 

is the Fisher matrix and the noise bias is 

b k = triC^QkC^N) . (7) 

Note that we could include the background power explicitly in these equations as a noise source 
when measuring the power in the Lya forest region, but we ignore this because its contribution is 
too small to change the weighting significantly. We will subtract it from the estimates later. The 
noise subtraction term, b k , is computed using the pipeline noise estimates, a p (not a w ), with the 
amplitude corrected as discussed above based on the differences between exposures. In principle, 
S in these equations should be the true covariance matrix; however, as we discuss below, we use 
the measured covariance from a previous iteration of the power spectrum determination instead. 

Except in a few cases that we will identify as they arise, when we set out to measure the power 
in a defined rest-wavelength region (e.g., 1041 < A rcs t < 1185 A for the Lya forest region) we first 
use equation (5) to estimate the power separately in halves of the region in each spectrum (e.g., 
1041 < A res t < 1113 A and 1113 < A rcst < 1185 A). Our choice of half-spectra is a compromise 
between competing desires for resolution in redshift and wavenumber. The full length of the forest 
in a spectrum corresponds to a redshift interval, Az ~ 0.4, that is unnecessarily large. While the 
precision of the measured power spectrum would support smaller than half-spectrum chunks to 
give finer redshift resolution than Az ~ 0.2, the shorter chunks would limit the /c-space resolution. 
Note that we could have used full chunks and still achieved the same z-resolution by more carefully 
applying the estimator equation, as we discuss below, but this would increase the computational 
time without much improvement in the final errors on the scales of relevance. 

After computing estimates P k for each half-spectrum, we perform a weighted average to de- 
termine Px(k, z) in redshift bins centered at Zi = 2 + 0.2 i where i = 1..13 (in this paper we only 
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present results up to i = 9). Each bin is the average of the power in all the half-spectra for which 
the redshift of the central pixel falls within ±0.1 of the bin center (we discuss below how we correct 
for an asymmetric distribution of data within a bin) . We combine sets of estimates using the Fisher 
matrices (equation 6) for the weighting. In practice this means that we sum the quantity in paren- 
theses in equation (5) over all estimates and multiply the result by the inverse of the sum of the 
Fisher matrices for each individual estimate. Our procedure would be optimal for Gaussian data, 
which the Lya forest is not; however, when we use the Gaussian approximation to compute the 
errors on the measured power the results are not much different from the more accurate bootstrap 
errors (see §3.1.2), so we conclude that our method is not far from optimal. 

Whenever we have a finite length of spectrum, there will be mixing between the power in 
different bins. Variable noise or gaps in the data will produce more mixing. This mixing is 
described in terms of a window matrix, which is given by the Fisher matrix in equation 6. In our 
standard procedure, the power spectrum estimates in equation 5 are multiplied with the inverse 
of Fisher matrix and are thus deconvolved with the window, which removes the mixing of other 
modes into the bin one is estimating (however, the bins are still correlated). This method thus 
produces a diagonal window matrix, so each combined estimate of Px(k,z) represents exactly the 
range of k corresponding to its bin. Our tests below show that there is no practical problem with 
instability in the Fisher matrix inversion (the window matrices are close to diagonal to begin with). 
A diagonal window matrix is desirable from a theoretical standpoint because our ability to compute 
the power spectrum from simulations is limited at both low k (by limited box size) and high k (by 
simulation resolution and complexity of physics). In the few cases where we use the power without 
deconvolution, we are using the estimator N k (F ~P x )k, where N k = (£) fc , F kW ) 1 (Seljak 1998). 

To compute the weight matrix C, we need an estimate of S, i.e., the power spectrum we are 
trying to measure. We solve this problem by computing Px(k,z) iteratively. The first estimate 
is made assuming S = 0. In subsequent iterations we compute S from the previous estimate of 
Px(k, z). This procedure converges quickly (the difference between S = and a reasonable estimate 
of S is significant, but once S is in the right ballpark it does not matter what it is exactly). We 
add a large constant (10.0) to all elements of the weight matrix, to remove all direct sensitivity of 
our power measurement to the mean of the chunk. This makes very little difference to the results 
on the scales we present. We are however still sensitive to the mean estimate from when we divided 
the spectrum by it. Even if the mean estimates are correct on average, the statistical error on the 
mean for each spectrum can still lead to a bias. If the errors on the mean estimate were small and 
uncorrelated with the fluctuations in the flux field, the bias would be 1 + 3 cr^, where a m is the 
error on the mean [to lowest order in a m , i.e., the bias is (l/(l + 5 m ) 2 ) ~ 1 + 3 (<5^), where 8 m is 
the fractional error in the mean, and (<5^) = We divide each estimate by this factor as part 
of our standard procedure; however, as we discuss below when we test our code on mock spectra, 
this approximation is not sufficient and we will need to include another small, fc-dependent, factor 
determined numerically using the mock spectra (this is the only use of the mock spectra other than 
for testing). 
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The reader may at this point be wondering what redshift the resulting Px(k,Zi) should be 
taken to represent, i.e., Zi is not necessarily the center of weight of the data, and neither is the 
mean redshift of the pixels in the bin, considering the rather complicated weighting in equation 
(5). In fact, the effective redshift is not even the same for each /c-bin in the same z-bin. We resolve 
this question - Px(k,z{) represents the power spectrum at precisely Zj (to first order) - in our 
construction of S ab = {S x 5 b x ) and = dS ab /dPk,i, where a and b label pixels at redshifts z a and 
Zb, and i labels the redshift bin in which this chunk of spectrum falls. To account for the evolution 
from z a and Zb to Zi, we define a power spectrum growth factor, D^^z) = [(1 + z)/(l + Zi)] ak '\ 
where 

dln[P x (k,z)} 
" M d\n{l + z) 



ln(P M+1 /n,*-i) (g) 



ln[(l + Zi+i)/(l + Zi_i)] 

(we use a one-sided derivative estimate instead of equation 8 for the first and last redshift bins). 
Now Qt b = Dki(z a b) Qt b \ , where z a b = (z a + z&)/2 and Qt b \ is computed as if the pixels were 
located at the center of the bin. Finally, S ab = ^2uQ ( £'Pk{ z i)- This correction may be difficult to 
understand intuitively at first, but it is really quite simple. The modification of Q just corrects 
the power spectrum estimate for the excess (dearth) of power that we expect for pixels in the high 
(low) redshift ends of the bin. The correction to S affects the weighting, simply producing a more 
accurate S at the redshift of the pixels in question. 

Note that an alternative method would be to treat the points Px(k,zi) as simply parameters 
of a continuous power spectrum defined by some form of interpolation. This would mean S ab would 
have non-zero derivative with respect to more than one of the power spectrum bins (e.g., usually two 
for linear interpolation). This method would be elegant, and probably produce narrower effective 
window functions in the z direction; however, it will increase the correlation in the z direction 
between measurement errors, because the same pixels would contribute to more than one power 
spectrum point. Since this more sophisticated method would allow long chunks of spectra to be 
used without degrading our z resolution, it would be most useful if we were trying to measure the 
power on even larger scales. 

How does our method compare to the straightforward Fourier transform (FT) method? The ba- 
sic FT method is to project the data vector, Sf, onto a set of modes of the form d a p = exp(i k a Avp), 
and to simply compute the variance of the amplitudes of all the modes with k in some bin, i.e., 



<\k a \<k 

! I 

P FT OC ^ 



^2 da P 5 f,p 

/3 



(9) 



where k m \ a and fc max define the bin, and the discrete spacing of k a is somewhat arbitrary (the 
natural spacing is Ak = 2ir/L, where L is the length of the spectrum, but nothing prevents one 
from choosing more finely spaced ks). Our estimator, equation (5), can be cast in a similar form, 
i.e., as a projection of the data vector onto a set of modes, and a sum of the squares of the mode 
amplitudes. We require that the mode amplitudes are statistically independent, which makes their 
computation equivalent to a computation of Karhunen-Loeve eigenmodes (see, e.g., Tegmark et al. 



-25- 



(1997)). Figure 7 shows the two most important modes for our bin with 0.00126 (kms~ ) _1 < k < 
0.00158 ( kms -1 ) -1 , for the chunk of spectrum shown in Figure 3a. In this case two modes differing 
primarily by a phase shift, analogous to sin(fcx) and cos(kx), contain most of the information, 
because our bin width is approximately 2ir/L. We see that the difference between our modes and a 
simple sine wave is not dramatic - there is a little bit of edge tapering (downweighting the edges to 
make the effective window on the data more compact in Fourier space) and some straightforward 
downweighting of the most noisy pixels. Curiously, there seems to be an additional effect where 
pixels adjacent to an edge are given extra weight, possibly as a way of compensating for missing 
data (this is seen more clearly in spectra where a narrow gap is present in the middle of the data). 
The picture is similar for bins with larger k, except of course that there are increasingly many 
important modes as the width of our bins increases (the bins have a fixed width Alog(fc), but the 
relevant mode width is Ak). For more discussion of the quadratic estimator see, e.g., Tegmark 



The method we adopt is optimal for Gaussian fields and therefore guarantees that no other 
method can surpass it. An additional advantage is that within this formalism window and covari- 
ance matrices are automatically computed. For continuous spectra with few gaps and near uniform 
noise one does not necessarily expect an FT method to be significantly worse. In practice the noise 
level is slowly varying across the spectrum, so averaging all the pixels uniformly is not optimal 
and degrades performance. Another advantage is that with our method each pixel pair has its own 
effective redshift and the correlations for a given pair are then interpolated to the redshift of inter- 
est using the appropriate evolution. In the FT method the whole spectrum is Fourier transformed 
first, so the redshift information is preserved only in an averaged sense, but a priori it is not clear 
how this average is defined. 



While the Fisher matrix obtained during the estimation process would give the error matrix 
for Px(k, z) if the data were Gaussian, we cannot reliably assume this. Our solution is to compute 
a bootstrap error matrix by the standard procedure (Press et al. 1992). From our data set of N 
spectra, we form a bootstrap data set by selecting N spectra at random, with replacement. The 



covariance matrix of Px(k,z) is taken to be = (AP x AP x \, where AP X = P x — (P x 



P l x is an estimate of the power in the zth bin from a bootstrap data set, and () means average 
over bootstrap realizations. We generally use 4000 realizations, after checking that this produces 
convergence in the result. We assume that the error correlations extend only one bin off diagonal 
in the z direction, because the spectrum of a single quasar practically never contributes to non- 
adjacent bins. 



(1997). 



3.1.2. Bootstrap Error Estimation 




We have no compelling reason to believe that this method of computing the error bars will 
give rigorously correct results. Considering the large number of off-diagonal elements that must be 
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Fig. 7. — Black solid and red dotted lines show the two primary modes onto which the data is 
effectively projected when we estimate the power in our bin centered on k = 0.00141 (kms -1 ) -1 , 
for the spectrum shown in Figure 3a. The horizontal axis scale in the figure is arbitrary. For 
comparison, the dashed line shows a simple sine wave with k = 0.00141 (kms -1 ) -1 . 
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estimated, one worry is that a particular linear combination of the bins may accidentally vary very 
little in our data set, so it will appear to have an unrealistically small error. Our tests on mock 
spectra (§3.2.2) show no sign of this problem. Still, to be conservative we apply one tweak to M after 
it is computed, in an attempt to inoculate it against the possibility. We perform a singular-value 
decomposition on M, which produces a set of independent vectors and their variances. We then 
compute the variances of the same vectors under the Gaussian approximation, using the Fisher 
matrix. If the bootstrap variance is smaller than the Gaussian variance we replace it with the 
Gaussian variance. Finally, we transform back to M. The tests on mock samples described below 
give us confidence that our procedure is reliable. 

3.2. Tests on Mock Data Sets 

We validate our procedure as implemented in code by applying it to mock data sets. Many 
iterations of these tests were required to produce results that show no serious problems in the error 
estimation or the power spectrum estimation itself. Testing the results on realistically created 
mock data is absolutely essential for measurements of such high precision. In §3.2.1 we describe 
our procedure for generating the mock spectra. We test our bootstrap error estimates in §3.2.2. 
Finally, we test the power spectrum estimation procedure for systematic errors in §3.2.3. 

3.2.1. Generating Mock Spectra 

We generate mock spectra by combining the auxiliary information we have for each observed 
spectrum (e.g., our continuum estimate, noise estimate, sky estimate, etc.) with a simplified version 
of the Bi et al. (1992) model for the exp(— r) field, which results in realistic looking spectra. 

For each observed quasar we start with the term we divide by before computing the power 
spectrum, A q C(\ xes \) (1 + 8c) S(X) (see equation 2). We multiply this by exp(— r) (generated as 
described below), smooth the result using the resolution from the observed spectrum, and sample 
the result onto the observed grid of pixels. This produces a noise free version of the flux we would 
observe coming from this quasar. We add flux from the sky as estimated for the observed spectrum, 
and transform the total flux to the number of photons that would be expected in each pixel. We 
generate a Poisson deviate with this mean, add the appropriate Gaussian read-noise for each pixel, 
transform back to the original flux units, and subtract the sky flux estimate to obtain an observed 
(noisy) quasar spectrum. The results of this procedure for each observed quasar are written into 
files in the same format as the observed spectra, so exactly the same code can be used to measure 
the power in the mock spectra. 

To generate the exp(— r) fields we use a simple model that is arranged to give roughly the 
correct power spectrum as a function of k and z, and the correct mean absorption as a function of 
redshift. For each observed spectrum, we start by generating a Gaussian random field, t^o, on a 
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very long, relatively finely spaced grid (65536 cells with width 7kms 1 , to be precise), with power 
spectrum 

w- '^y ^"' <io> 

where ko = 0.001 (kms" 1 )" 1 , v = 0.7, and Rs = 5kms _1 [this Ps was chosen after some ex- 
perimentation because it produces a final flux power spectrum with approximately the same k 
dependence as the observed Pp{k, z)\. An arbitrary cell in this grid is chosen to correspond to the 
redshift of the quasar, and the evolution of the amplitude of the power spectrum with redshift is 
imposed by the transformation <5j = a(zi) <5j 5 o with a 2 (zi) = 58.6 [(1 + Zj)/4]~ 2 ' 82 , where the form of 
a(z) was chosen so that the final flux power spectrum would evolve like the observed one. Next we 
make the squared lognormal transformation n$ = [exp(<5j — cr?/2)] 2 , where af is computed from the 
input power spectrum, including the amplitude factor (the factor af/2 in the exponential just fixes 
the mean of the lognormal field to 1). We smooth the n field with a Gaussian filter with rms width 
R T = 20kms _1 and multiply it by a factor 0.374 [(1 + Zj)/4] 5 ' 10 to produce a field r (this redshift 
evolution factor produces roughly the observed redshift evolution of F). The mock transmitted 
flux in each grid cell is then Fi = exp(— n), which is sampled as described above. 

The procedure described above leads to realistic looking spectra of the Lya forest. We have 
verified that it generates a bispectrum that is within a factor of 2 of the one measured in N- 
body simulations. The main advantages of this procedure over the N-body simulation approach 
when generating the mock spectra are that it is faster, so one can make an arbitrary number of 
independent realizations, and that the simulated spectra can be of arbitrary length, important to 
eliminate any periodicity effects (this would be impossible with simulations, where a typical box 
size is much shorter than the total length of a single spectrum). Both of these advantages are 
critical for a high precision test. We determine the true Pp(k,z) by a simple FFT of extremely 
long exp(— r) fields (without redshift evolution). 



3.2.2. Tests of the Error Estimates 

Without accurate statistical errors it is difficult to identify systematic problems, so we first test 
our bootstrap procedure for estimating the errors. Note that there is no reason to expect bootstrap 
errors to be perfect (there is even some ambiguity in how exactly the bootstrapping should be done 
when the data do not consist of statistically identical objects). Regardless of systematic errors in 
the method, the only difference between the power spectra measured from two mock data sets that 
differ only in the random seed that was used to create them should be the statistical errors that we 
estimate. We test our error bars by generating ten different sets of mock data and computing \ 2 
for the differences between each of them and their error weighted mean, using the bootstrap error 
bars and the 108 points in 0.0013 < k < 0.02 (kms -1 ) -1 , and 2.1 < z < 3.9. This is effectively a fit 
of 108 parameters to 1080 data points, with 972 degrees of freedom. The total \ 2 ls 939, perfectly 
consistent with a random fluctuation around the mean, and strongly disfavoring an underestimation 
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of the errors by more than a couple percent. 



3. 2. 3. Tests of the Power Estimates 



We can now search for systematic errors. To enhance the statistical significance of any errors, 
we average our ten sets of mock spectra to form a single, more precise measurement. The result 
is shown in Figure 8. The results look reasonably good; however, we find an unacceptably bad 
X 2 = 346 for the comparison between our measured Pp(k,z) and the true power spectrum (there 
are 108 degrees of freedom). To quantify the systematic problem, we first assume the bias has the 
form B(k) = Bq [fc/0.0067 (kms -1 )" 1 ]" = P mea sured/ Pinput and fit for the values of Bq and v that 
minimize x 2 in the comparison. We find B = 1.0036 ± 0.0014 and v = -0.0173 ± 0.0013 with 
X 2 = 173 for 106 degrees of freedom [the pivot point fco = 0.0067 (kms -1 ) -1 was chosen to make 
the errors independent; the amplitude coefficient would be larger if we were not already dividing 
by 1 + 3 as explained in §3.1.1]. The combination of slope and amplitude errors corresponds 
to a 3.1% excess of power at our largest scale, k = 0.0014 (kms -1 ) -1 , and a 1.3% underestimate 
at k = 0.018 (kms -1 )" 1 . We find some less significant dependencies by generalizing the fitting 
formula even more to 



where a(z) = (1 + z )/(l + z), with z = 2.85. The parameters are B = 1.0073 ± 0.0016, 
H = 0.049 ± 0.012, v = -0.0195 ± 0.0015, rj = -0.0157 ± 0.0038, and C = -0.026 ± 0.012, with 
X 2 = 135. Where does this bias come from? We expect some bias related to the division of each 
chunk of spectrum by its overall mean (not because of an integral constraint suppression of large- 
scale power - our estimator should take care of that - but because of statistical error in the estimate 
of the mean that we divide by). When we measure the power without this division by the mean, 
which we can only do using mock spectra, we find significantly smaller corrections - small enough 
to ignore when model fitting. 

We expect that this bias should be present when we use real observed spectra, so we will 
correct for it by dividing the measured power by B(k,z). We describe its effect on the amplitude 
and slope of the power spectrum below (table 1). 



Figure 9 shows the raw power measured in our standard Lya forest rest wavelength range, 
1041 < A rC st < 1185. All the figures in this subsection show Pio4i,ii85i n °t the background sub- 
tracted power Pp. Our normalization convention is: 




(11) 



3.3. 



Raw Power 




(12) 
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Fig. 8. — Error bands show the average power spectrum [A 2 (A;) = ir^ 1 k P(k)] measured from 
ten sets of mock spectra. Lines show the true power. Redshift bins, strictly from bottom to top, 
are: black, solid line, z=2.2, blue, dotted, 2.4, cyan, dashed, 2.6, green, long-dashed, 2.8, magenta, 
dot-dashed, 3.0, red, dot-long-dashed, 3.2, black, thin, solid line, 3.4, blue, dotted, thin, 3.6, cyan, 
dashed, thin, z=3.8. 
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Fig. 9. — Error bars show the power spectrum measured from the observed spectra in the wavelength 
range 1041 < A ros t < 1185. The lines connect the points to identify them and to guide the eye. 
Redshift bins, from bottom to top (roughly) are: black, solid line, z=2.2, blue, dotted, 2.4, cyan, 
dashed, 2.6, green, long-dashed, 2.8, magenta, dot-dashed, 3.0, red, dot-long-dashed, 3.2, black, 
thin, solid line, 3.4, blue, dotted, thin, 3.6, cyan, dashed, thin, z=3.8. 
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We usually plot the dimensionless quantity A 2 (k) = tt 1 k P(k), the contribution to the variance 
per unit In k. 

Figure 10 shows the fractional errors on all of the measured points. The errors are less than 5% 
for most of the points, and frequently as small as 3%. If we were only estimating a single amplitude 
parameter by combining all these points then its error would be 0.6%. An overall logarithmic 
slope would have an error ±0.005. The errors on the largest scales are increased somewhat by the 
diagonalization of the window matrix. 

Figure 11 shows the ratio of subtracted noise power to measured signal power (-Pio4i,ii8s) f° r 
each point. The noise power is significant (20-30%) on all scales, but diverges at high k where the 
resolution suppresses the absorption power. The lowest redshift bin has the most noise, due to the 
lower Lya forest power combined with extra noise at the short wavelength end of the spectra. 

Figure 12 shows our window matrix (at z = 2.6), which we proceed to diagonalize. The matrix 
is reasonably close to diagonal already, with large contributions only from adjacent bins. It is useful 
to diagonalize the matrix at this stage, rather than waiting until the model-fitting stage, because 
this allows us to compute bootstrap errors directly for the final bins (the bootstrap error calculation 
and window matrix diagonalization do not perfectly commute). 

Figure 13 shows the ratio of the bootstrap errors to the errors estimated assuming the data 
are Gaussian. We did not apply the Gaussian floor to the bootstrap errors when making this 
figure. Typically the bootstrap errors are 0-20% larger than the Gaussian errors. Figure 14 shows 
examples of the estimated correlation between the errors, at z = 2.6. We note that diagonalizing 
the window matrix noticeably reduces the error correlations. The bootstrap errors are, in contrast 
to the Gaussian errors, noticeably correlated ((APjAP,) jop^opj ~ 0.0—0.2 when > 1, where 
i and j label the bins) across the full k range. These differences between bootstrap and Gaussian 
errors are not necessarily an indication of intrinsic non-Gaussianity in the absorption fluctuations. 
Possible alternative explanations for the differences include the uncertainty in the mean flux value 
that each chunk of spectrum is divided by and the uncertainty in the noise-subtraction term for each 
chunk, neither of which are included in the Gaussian estimate and both of which would increase 
the error in a way that is correlated across k bins. 



3.4. Background Subtraction 

Our background subtraction is the power in the wavelength range 1268 < A rcst < 1380 A. Fi gure 
15 shows -Pi268,i380 and fio4i,H85 for comparison. The bump at k ~ 0.013 (kms -1 ) -1 in Pi268,i380 
is probably due to the CIV doublet at separation 499 kms -1 . The bump at k ~ 0.003 (kms -1 ) -1 
may be due to the SilV doublet at separation 1933 kms -1 . Figure 16 shows fi268,i38o/-Pi04i,ii85- 
We see that, even though the background power is a small fraction of the Lya forest power, it 
is quite significant when compared to the small size of the errors on the Lya forest power. It is 
important to remember that even a small overall systematic error can be very significant if it covers 
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Fig. 10. — Lines connect the fractional errors on each measured Pio4i,il85(&> z ) point, using the 
usual line-type and color and scheme (see Fig. 9 - the highest two curves are the highest two 
redshifts, the lowest is z = 2.8). 
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Fig. 11. — Lines, with types following the usual pattern (see Fig. 9), connect the quantity 
-fnoise/-Pi04i,ii85 f° r each measured point (the highest line is the lowest redshift). 
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Fig. 12. — The window matrix for bands indicated by the maximum (before diagonalization) . 
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Fig. 13. — The ratio of bootstrap errors to the Gaussian estimate of the errors. See Fig. 9 for 
correspondence between lines and redshift bins. 
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Fig. 14. — Examples of the correlations between the errors, (APiAPj) /crp t i crpj. The black solid 
lines and squares show the error correlation when the window matrix is diagonalized. The red 
dotted lines and triangles show the correlations between points before diagonalization. The lines 
marked by symbols are the bootstrap estimate, while the unmarked lines are the Gaussian estimate. 
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0.001 0.01 



k [(krn/s)- 1 ] 

Fig. 15. — The upper set of lines show Pio4i,ii85> the lower set of lines show Pi268,i380- The colors 
and line types identify redshift bins as defined in Fig. 9. 
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k [(krn/s)- 1 ] 

Fig. 16. — The lines show the ratio -Pi268,i38o/-Pi04i,ii85- The uppermost (black solid) line is z = 2.2, 
and the next (blue dotted) is z = 2.4 (see Fig. 9 for the rest of the line definitions). For reference, 
the error bars starting at zero show the fractional errors on Pio4i,H85(&> z = 2.6), which are much 
larger than the errors on Pi268,i380 (we are simply plotting ap(k)/P(k) as in Fig. 10, except that, 
for clarity, we show error bars starting at zero instead of a connected line). 
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many data points (e.g., a 1/2 a error over 100 points shifts the mean by 5 a). 

We are going to subtract the power in the range 1268-1380 A from the Lya forest power, but 
it is informative to measure the power at other places in the quasar rest frame for comparison. 
The range 1409-1523 A includes CIV absorption (at 1548.2 and 1550.78 A) but excludes SilV (at 
1393.75 and 1402.77 A) and shorter wavelength transitions. Figure 17 shows -Pi409,i523/-Pi04i,ii85- If 
all of the power was coming from metal line absorption, the power in the range 1409 < A res t < 1523 
A should always be less than the power in the range 1268 < A rcst < 1380A. As we see in Figure 
18, which shows the difference in the background fractions, (-Pi268,i380 — -P1409, 1523 )/-Pio4i,i 185 j the 
power in Pi268,i380 is greater than ^1409,1523 except on large scales. The difference on large scales 
suggests that there is tiny amount of power left in the quasar continua (in spite of our division 
by the mean continuum), which is larger in the range 1409-1523 A than in the range 1268-1380 
A. Finally, Figure 19 shows i 3 i558,i774/-F > i04i,ii85> P a st the wavelength of CIV absorption. The 
reduction of power relative to shorter wavelengths is dramatic, but not surprising since CIV is the 
most common metal absorber. It does suggest however that most of the power is due to metals 
and not continuum fluctuations, unless the continuum in the range 1558-1774 A has significantly 
less power relative to other intervals studied here (which is admittedly not inconceivable). It 
seems likely, although we are not certain, that the z = 2.2 background power has a noticeable 
contribution from measurement-related problems, because the alternative is a very sudden increase 
in metal absorption power. 

What is the upshot from these studies? The metal absorption appears to contribute a small 
but significant amount of power, which should also appear in the Lya forest region. We subtract 
this power from the power measured in the forest. There is some indication of measurement-related 
problems in our lowest redshift bin. The power contributed by deviations of the quasar continua 
from their mean appears to be small. 

While the idea that -Pi268,i380 contains almost exclusively power due to simple metal absorption 
seems plausible at this point, when we perform consistency checks in §4.4 we find evidence that 
this is not the case. Splitting the data set used to measure -Pi268,i380 i n half based on the noise 
level in each spectrum, we find that the power in the halves is significantly different, by as much 
as 50% in some bins. Splits based on several other properties of the data (e.g., sky to quasar flux 
ratio) also show significant differences, but we find that these differences can all be accounted for 
by the difference in the basic noise level in the subsamples. Splits of the Lya forest data set show 
similar trends in Pio4i,ii85 with the splitting parameters, although the fractional differences are 
much smaller. While we don't know the source of this noise dependence, it is not hard to imagine 
relatively benign reasons for it. For example, if sky subtraction is imperfect this would add an 
increasing amount of power as the sky flux, and thus noise level, increases relative to the quasar 
flux. The procedure we describe next would remove this power. 

Since we know that -Pi268,i380 depends on noise it seems logical to subtract the value of -Pi268,i380 
corresponding to the level of noise in the forest, rather than the best measured value of -Pi268,i380> 
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Fig. 17 — The ratio of Pi 4 09,i523 to Pi 4i,ll85- 
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Fig. 18. — The difference between -Pi268,i380 and Pi409,i523, divided by Pio4i,ii85, with the fractional 
errors on Pio4i,ii85(&> z = 2.6) plotted as usual. 
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Fig. 19. — The ratio of -Pi558,i774 to Pio4i,ii85- The uppermost (black solid) line is z = 2.2. 
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which is dominated in practice by data with less noise. If we had simply misestimated the noise 
by an overall factor, for example, the errors in Pio4i,ii85 an d -Pi268,i380 would cancel for this form 
of subtraction. To implement this idea, we model the background subtraction term as a linear 
function of the noise level, 

^I268,i38o0, z, a w ) = A(k, z) + B(k, z) a w , (13) 

where a w is the mean noise level in the data computed in the same way as the mean flux level (this 
is the mean of the normalized noise, i.e., after division by continuum and mean flux). The choice 
of a linear relation is arbitrary but it does the job (see §4.4) better than the alternatives we tried. 
We compute A(k,z) and B(k,z) for each value of k and z using a linear fit to the full sample of 
spectra that probe -Pi268,i38o(fc> z), weighted by the Gaussian estimate of the error on each point for 
each spectrum. When the time comes to subtract the background from Pio4i,ii85 to obtain Pp, we 
use a w computed in the 1041 — 1185 A wavelength range to compute the appropriate subtraction 
term. Figure 20 shows the extra power subtracted through Equation (13), beyond what we would 
subtract if we simply used -Pi268,i38o(&> z) from Figure 16. It is typically less than 4% of the Lya 
forest power, but rises to 10% at the highest k for the lowest redshift. 

The reader who is paying attention may complain that we have no compelling reason to believe 
that this source of noise-dependent power that we do not understand depends on noise in the same 
way inside and outside the Lya forest region. This would be true, except that when we follow 
this prescription for background subtraction the differences in Pio4i,ii85 between subsamples are 
removed (see §4.4). This would be a remarkable coincidence if our model for the subtraction was 
not substantially correct. 

Note that the background power has much smaller (absolute) statistical errors than Pio4i,ii85) 
mostly because there is simply less power, but also because there are more quasars probing a given 
redshift interval. 

4. Consistency Checks 

In §4.1 we describe how we use fits of theoretical models to the Pp(k,z) results to help un- 
derstand the importance of any systematic errors. We plot the correlation function of the Lya 
forest in §4.2 and use it to identify a significant contribution to Pp(k,z) from Silll absorption. In 
§4.3 we examine the effects of modifications of our procedure. In §4.4 we break the data set up in 
many different ways to look for dependencies of Pp(k,z) that should not exist. In §4.5 we discuss 
the possibility that continuum fluctuations contribute significant power. Finally, we compare our 
results to past measurements in §4.6. 
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Fig. 20. — The difference between the noise dependent background power that we subtract through 
Equation (13) and -Pi268,i38Cb relative to -Pio4i,ii85> i-e., this is the extra fractional power that is 
subtracted because we correct for the difference between the typical noise level in the forest and in 
the range 1268-1380 A. 
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4.1. Rudimentary Fitting of Theoretical Models 

The ultimate purpose of measuring the Lya forest power spectrum is to determine cosmological 
parameters by comparing the observed Pp(k, z) to the predictions for different cosmological models. 
For the ACDM models supported by present observations, the universe is nearly Einstein-de Sitter 
at z > 2, so cosmology influences the Lya forest almost entirely through the linear theory power 
spectrum of the mass fluctuations, Pl(A;, z), at z ~ 3 and k ~ 0.01 (kms" 1 )" 1 (roughly 1 comoving 
h/Mpc, depending somewhat on the model). We usually parameterize Pi(k,z) by its amplitude, 
A 2 (k p ,z p ) = kpPL(k p ,z p )/2iT 2 , slope n e g{k p ,z p ) = (ilnP^/dln k\ kp , and curvature a e s(k p ,z p ) = 
dn e ff/dlnk\ kp , where we use k p = 0.009 (kms" 1 )" 1 and z p = 2.6 as the pivot points. 

A full discussion of the details of theoretical modeling of Pp(k, z) using numerical simulations 
is beyond the scope of this paper. Furthermore, the theory of the Lya forest is perhaps less 
certain than the observations, so we want to present the observational results un-tarnished by 
theoretical interpretation. However, it is very useful to interpret the possible systematic errors in 
the appropriate context of cosmological model fitting. In other words: without model fitting, it is 
difficult to know how important a given change in Pp(A;, z) is. 

In this paper we take a cautious approach to the theoretical model fitting - we perform fits to 
different estimates of Pp(k, z) computed using modifications of the extraction procedure or different 
subsamples of the data, however, we do not give the central result, only the deviations in the results 
from the value obtained from our preferred Pf(/j, z). These deviations in fitting results should give 
the reader a useful indication of the importance of systematic effects in the data, regardless of the 
reader's opinion of the theory. 

The simulations and fitting procedure that we use are described in McDonald et al. (2004), 
where we present the final result. We use a ACDM transfer function, and perform the fit with A 2 
and n e fr as free parameters (because a e ff = dn e g jd In k is not tightly constrained by the present Lya 
forest data alone, we fix the primordial running a = dn/dlnk, not to be confused with a e g ~ —0.2, 
to zero). Unless otherwise specified, we perform the fits using the 108 Pp(fe, z) points in the ranges 
0.0013 < k < 0.02 (kms -1 ) -1 and 2.1 < z < 3.9. We allow for some error in our noise estimate 
by permitting the noise subtraction terms to vary independently in each redshift bin by 5% (9 
extra free parameters to fit for, constrained by Gaussian likelihood function with this rms width). 
We also allow a single overall parameter describing the squared resolution error to vary with rms 
constraint (7 kms" 1 ) 2 . 

The Lya forest model in the simulations is controlled by the externally constrained functions 
F(z), the mean absorption, T±^(z), the temperature at overdensity 1.4, (7 — l)(z), the logarithmic 
slope of the temperature-density relation, and a reionization parameter that we will call x re i. F{z) 
is described in our fits by the 10 parameter formula Fj = TFi, where i labels our 9 redshift bins, Fj 
gives the arbitrarily normalized z dependence and T is an overall normalization. We have performed 
a preliminary analysis using the formalism in §2.5 to measure F(z) from SDSS data and we use 
this to constrain the parameters Fj (the error on each redshift bin is ~ 0.005). Because the SDSS 
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analysis does not constrain the overall normalization, we leave T free except for the additional 
constraint that we require Fi interpolated to z = (3.9, 3.0, 2.4) to match the HIRES constraints 
F = (0.458 ± 0.034, 0.676 ± 0.032, 0.816 ± 0.023) (see McDonald et al. (2000) - we have modified 
the numbers slightly and increased the errors to allow for systematic uncertainties, as discussed 
in Seljak et al. (2003)). We parameterize T L4 (z) and (7 — l)(z) by quadratic functions of z (3 
parameters each) with the external constraints T 1A = (20100 ± 3400, 20300 ± 2400, 20700 ± 2800)K 
and 7 - 1 = (0.43 ± 0.45, 0.29 ± 0.3, 0.52 ± 0.14) at z = (3.9, 3.0, 2.4) (see McDonald et al. (2001) 
- we added 2000 K in quadrature to the temperature errors to allow for systematic errors). Note 
that there are other, sometimes more precise, measurements of F (Schaye et al. 2003; Bernardi 
et al. 2003) and the temperature-density relation (Schaye et al. 2000; Ricotti et al. 2000) in the 
literature - our choice of McDonald et al. (2000) and McDonald et al. (2001) for this example is 
simply for convenience. The redshift of reionization and post-ionization temperature of the gas 
influence Lya forest predictions because the smoothing of the gas on small scales depends on its 
pressure history. At the level of precision we care about, this dependence can be captured by a 
single parameter. In our modeling, we use x re i to interpolate between two reasonable boundaries, 
reionization heating of the gas to 25,000 K at z = 7 or to 50,000 K at z = 17, both of which are 
consistent with our temperature constraints T\^{z). However, in this paper we fix x re «, because it 
is weakly constrained by the data and the hard lower limit we have to impose on the redshift of 
reionization leads to non-Gaussian errors on the power spectrum parameters we are interested in 
(this is a problem of presentation, not of principle). 

Figure 21 shows our first fit to the standard Pp(k,z) results. The value of \ 2 = 193.7 is 
much too high for approximately 106 degrees of freedom (we are marginalizing over a large number 
of nuisance parameters, but these generally are externally constrained so they do not reduce the 
number of degrees of freedom). Including a e s as a free parameter does not improve the fit signifi- 
cantly. It appears that much of the disagreement comes from bumps in the power spectrum, e.g., 
at k ~ 0.003 (kms^ 1 ) -1 . This motivates us to look at the correlation function of the flux. 



4.2. The Correlation Function and the Silll Cross-Correlation 

Sometimes features in the power spectrum are easier to understand by looking at the correlation 
function, £(v) = (5(x)5(x + v)) (v is as usual a stand-in for wavelength differences, as is x in this 
case). We show the normalized correlation function, £(v)/£(0) for the first six redshift bins in Figure 
22. The correlation function shows the expected behavior - positive for small v, negative for large v 
- except for an obvious bump at v ~ 2200 km s" 1 . We focus on this bump in the inset panel of Figure 
22. The most likely explanation seems to be cross-correlation between Lya and Silll absorption. 
Silll absorbs at A = 1206. 50A, so the Silll absorption by gas at some point along the line of sight 
will appear in the spectrum separated by 2271 km/s from the Lya absorption by the same gas. 
We see that the bump in £(v) appears at this separation, and note that the features that ruined 
our power spectrum fit appear at the expected multiples of k = 2vr/2271 = 0.0028 (kms -1 )- 1 . 
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Fig. 21. — Points with error bars show the observed Pf(/s, z). Lines show our first attempt to fit a 
theoretical model, which is not a good fit to the data. From bottom to top — z=2.2: black, solid 
line, open square; z=2.4: blue, dotted line, 4-point star (cross); z=2.6: cyan, dashed line, filled 
square; z=2.8: green, long-dashed line, open triangle; z=3.0: magenta, dot-dashed line, 3-point 
star; z=3.2: red, dot-long-dashed line, filled triangle; z=3.4: black, thin solid line, open pentagon; 
z=3.6: blue, thin dotted line, 5-point star; z=3.8: cyan, thin dashed line, filled pentagon. 
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Fig. 22. — The normalized correlation function, £(v)/£(0), in the Lya forest region, uncorrected 
for resolution. In the inset panel, the solid lines show an expanded view of the Silll-Lya cross- 
correlation bump, the dashed line shows 0.04 £(v — 2271)/£(0) for comparison, and the vertical 
dotted line marks v = 2271 km s" 1 . Note that there is no evidence for any other metal with 
wavelength close to Lya transition being important. In particular, we see no bump at ~ 5600 kms -1 
or ~ 6700 km s -1 , corresponding to NV-Lya velocity differences. This correlation function should 
not be used for any quantitative science, as we have not corrected for resolution effects, have not 
checked carefully for systematic errors, and have not given statistical errors. 
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Note that this is the only correlation seen; another metal correlation one might expect to see is NV 
(A = 1238.8, 1242. 8A), but there is no apparent feature at the corresponding velocity differences 
(Av ~ 5600, 6700km/s), as seen in figure 22. 

What should we do about this Silll-Lya cross-correlation, since the poor x 2 suggests that it is 
too significant to ignore? Our first guess might be that the Silll-Lya correlation is a simple offset 
version of the Lya-Lya correlation, i.e., something like Csiiii-Lya('f) oc £Lya-Lya(M — 2271kms~ 1 ). 
The simplest way to model this is to assume that the Silll structure is equal to that of the Lya 
forest up to an overall normalization, 5f = S(v) + a 5(v + v 3 ), where S(v) is for Lya only and 
v% = 2271 km s _1 . The corresponding correlation function is 

£ F (v) = (1 + a 2 ) +a^(v + v 3 ) + a £(w - v 3 ), (14) 

with corresponding power spectrum 

P F (k) = (1 + a 2 ) P(k) + 2 a cos(v 3 k) P (k) , (15) 

where unsubscripted £ and P are understood to mean Lya-Lya. For our first fit to Pp{k,z) 
accounting for Silll using equation (15), we assume a = f/[l — F(z)], with / as a single extra 
free parameter of the fit. We find a remarkable improvement in x 2 , from 193.7 to 130.9. We 
find / ~ 0.011 (a ~ 0.04, depending on the redshift). The small value suggests that the relative 
contribution of Silll to the autocorrelation is a 2 < 0.004, which will not affect our fit results 
significantly (see §4.3). We thus only need to estimate the cross-correlation term. We also tried 
allowing a power law 1 + z dependence for /, but the improvement in fit, A% 2 = 1.1, was not 
significant. 

In the inset panel of Figure 22 we plot scaled £(v — 2271), to show how the shape of the bump 
compares to the zero-lag correlation. It is difficult to compare the shapes by eye, because of the 
slope of the underlying correlation, but it appears that this model explains the cross-correlation 
reasonably well. We can allow for a change in scale using the slightly more general form 

Zf(v) = + a£[s(v + v 3 )] + a £ [s (v - v 3 )} . (16) 

Allowing s to vary freely only improves x 2 by 0.7 (note that the logarithmic fc-binning that we 
use may reduce our ability to constrain these parameters). The error bars on other parameters 
may increase when we include z dependence of / and allow s to be free, so to be conservative one 
probably wants to leave them free even though they are not needed. In our standard fit in this 
paper, we allow / to have z dependence but fix s = 1. We show the improved fit to Pp{k,z) in 
Figure 23. 

4.3. Modifications of the Procedure 

The pipeline developed for this analysis includes many improvements and corrections that were 
added throughout the development. It is worth taking a step back to ask how important the various 
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Fig. 23. — Points with error bars show the observed Pp(k, z). Lines show our best fit after including 
Silll absorption approximately in the theory. See Fig. 21 for line and point definitions. 
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corrections are for the final result. Table 1 lists various modifications of our procedure (described 
individually below), and quantifies their effects on the fit results. In each case we show the change 
in the best fit A 2 and n e s relative to our standard fit, and their error bars for comparison to the 
standard fit. We give x 2 to indicate the goodness of fit of the theory to the modified measured 
Pp{k,z). We reiterate that we are not asserting the correctness of the theory that we use in the 
fitting - we give these \ 2 values and other fitting results only to show trends. We list A% 2 between 
the standard procedure best fit and the variant best fit, to give an indication of how significant the 
deviation is (this is necessary because the errors are correlated so simply knowing A 2 and n e s and 
their errors does not give the full picture - see Figure 26 for an example of the full error contours). 
Because the statistical fluctuations between these different fits should be small, a 1 a difference 
(or even less) should be interpreted as "significant", however, since we believe that our standard 
fit is better or more conservative than all of the variants, our systematic error should generally 
be smaller than the deviation shown. Note that where applicable the changes in procedure are 
only applied to the fio4i,H85 calculation, not the -Pi268,i380 result that is used in the background 
subtraction (small changes in -Pi268,i380 have no effect on the final results). 

Our first variant is to not diagonalize the window matrix. Figure 24 shows the measured power 
spectrum before and after diagonalization. Figure 25 shows the ratio of the diagonal errors after 
diagonalization of the window matrix to before diagonalization. Not diagonalizing the window 
matrix does lead to a significant change in the fitted parameter values, and the error on n e s 
decreases by about 12%. We are, in effect, using information from a wider range of scales, but this 
forces us to use theory results outside their range of validity (e.g., at low k we need to extrapolate 
beyond the size of our simulation boxes). Note that the change in error on n e g, from 0.024 to 
0.021 implies that we should expect a random difference between the two results with typical size 
(0.024 2 — 0.021 2 ) 1 / 2 = 0.012, i.e., what might seem like a surprisingly large part of the difference 
between the results could be random. 

As discussed above, the correction for Silll-Lya correlation is very important to the goodness 
of our fit. It is less important for the best fit values, changing them only by 0.8 a for n e ff and 0.4 a 
for A 2 . Normally we allow a power law dependence on redshift for the amplitude of the Silll-Lya 
correlation, but removing this freedom makes almost no difference. Allowing the correlation scale 
for the Silll-Lya correlation to be different than for Lya-Lya (freeing s in eqn. 16 - we usually 
fix this in this paper for technical reasons) has only a very small effect. Including the Silll-Silll 
autocorrelation term (the a 2 part of eqn. 15) in the fit has essentially no effect. 

For our standard fit, we allow variation in the noise amplitude at each redshift, represented by 
a multiplicative parameter subject to a 5% rms Gaussian constraint. Our fitting procedure then 
marginalizes over this component. Reducing this constraint to 0.5% (i.e., fixing the noise) produces 
no change in our fit result, and does not even reduce the error bars noticeably. Leaving the noise 
essentially free makes a noticeable difference in the fit results, decreasing the amplitude by about 
2/3 <7, increasing its error by 20%, and decreasing x 2 to 123.8. Changes at this level are expected 
when we remove the constraints on some parameters, and do not imply that the constraints were 
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Fig. 24. — Dotted lines connect the power spectrum points before diagonalization of the window 
matrix. Solid lines show the points after diagonalization. 
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Fig. 25. — Ratio of the diagonal errors on our Pio4i,ii85 {k, z) points (connected) after the window 
matrix diagonalization to before the window matrix diagonalization, with the usual line definitions. 
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too small (i.e., we are effectively removing 9 data points from the fit so we generally expect a 
decrease in x 2 ■, increase in the error bars, and some corresponding drift in the parameter values). 
Removing our spectrum-by-spectrum noise estimation makes very little difference. Finally, we note 
that if we did not correct the noise as discussed in §2.3, the results would change significantly. As 
an example, we show the fit results in the 6 In A 2 — 5n c g plane in Figure 26, for our standard case 
and these noise-related variants. We show the ratio of the power without individual noise estimates 
for each quasar to our standard case in Figure 27. 

Our requirement that the principal components of the error matrix have at least the Gaussian 
level of variance makes no difference to the fit results, although it improves x 2 f° r the fit a little bit. 
Simply using the Gaussian error matrix instead of bootstrap errors makes no difference to the fit 
results but increases % 2 significantly. Ignoring the bootstrap error correlations increases the error 
on n e s by about 12%, and significantly reduces x 2 ■ 

We normally use the range 1268 < A rcst < 1380A for our background subtraction (i.e., subtract 
Pi268,i38o)- Removing the background subtraction entirely reduces the inferred amplitude by 2a, 
and the slope by la, and results in a very large x 2 (the error on n c g- also decreases significantly, but 
this is mostly because of the change in the best fit values, not because of uncertainty in the back- 
ground subtraction). (Note that removing the background subtraction, which increases Pp{k,z), 
decreases the inferred amplitude because the fitted F decreases more than enough to offset the 
increase in Pp(/s, z).) Clearly the background cannot be ignored. Using -Pi409,i523 instead produces 
a somewhat disturbingly large 0.028 (1.1 a) increase in n e ff. We expect the longer wavelength range 
to give a less accurate estimate of the background power, because some metals are missing, but 
further investigation shows that most of this difference is probably caused by CIV BALs adding 
power to the 1409 < A rcst < 1523A region. As we see in Table 1, removing the adjustment for noise 
dependence of the background (see Equation 13) brings the two background regions closer together 
(this is reflected in Figure 18). Adding the 147 BAL quasars identified by our automated algorithm 
leads to a huge discrepancy (0.094 in n e g) when we use the -Pi409,i523 background, but only when we 
adjust for noise level (without this the discrepancy for n c g, not shown in the table, is only 0.029). 
Note that the BAL cut makes essentially no difference to our standard fit using the -Pi268,i380 back- 
ground. All of these differences are easy to understand: First, BALs are known to be strongest in 
CIV absorption (Hall et al. (2002)), so it is not surprising that we see the effects of BALs primarily 
in this wavelength region. Second, both our original by-eye and subsequent automated removal 
of BALs inevitably identify the features more easily in less noisy data, so the power from BALs 
naturally shows up when we intentionally use the noisier spectra for the background power. The 
fact that removing the 147 most obvious BALs has essentially no effect on our basic result gives us 
confidence that any remaining BAL features in the Lya forest and 1268 < A rcst < 1380A regions 
are not significant. 

To investigate the effect of a systematic uncertainty in the spectral resolution, we include in 
our fits an overall factor of the form exp(afc 2 ) multiplying the power spectrum, where a is a free 
parameter. In our standard fit we impose an external constraint on a, ±(7kms -1 ) 2 rms. Essentially 
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Fig. 26. — Fit results for variant noise treatments. The error bars show the 1 a error on each 
parameter. The ovals show A% 2 = 2.3. Standard case (with 5% noise amplitude freedom in each 
redshift bin): black, solid lines. No individual noise estimate for each quasar: red, dotted lines. 
Noise amplitude freedom 0.5% (50%): green, short-dashed lines (blue, long-dashed lines). The 
magenta, dot-dashed line shows the result using the pipeline noise estimates. 
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Fig. 27. — Ratio of Pio4i,ii85 computed without quasar-by-quasar noise re-estimates (a constant 
16% extra noise power was assumed instead) to the standard case. 
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removing this freedom has no effect on the fit, while leaving a essentially free increases the error 
on the amplitude by 40%, and increases n e g by 2 /3a (the change in fitted amplitude is certainly 
consistent with drift from the increased error). As we show in Figure 5, our standard fit should be 
conservative. 

Simply dividing each chunk of spectrum by its mean instead of also dividing by the mean 
continuum before estimating the power from the chunk makes little difference to the fit results. 
Division by the mean continuum actually increases the measured flux power by ~ — 2%, as we 
show in Fig. 28. We have performed a preliminary PCA analysis to try to model fluctuations around 
the mean continuum. When we use continua for each quasar composed of 13 PCA eigenvectors, 
our results change only a little (n e g is reduced by 0.4<r), and x 2 increases, probably an indication 
of the unsatisfactory level of noise that we know remains in our estimates. As we see in Figure 29, 
the modification of adding a large constant to the weight matrix to make our measurement less 
sensitive to the mean of each chunk has little effect (the effect is larger on larger scales). 

The line "no bin-redshift correction" in Table 1 refers to removing the correction for evolution 
in the power across the width of the redshift bins (see eqn. 8). We see (Fig. 30) that this correction 
mostly affects the lowest redshift bin (where the \ow-z edge of the bin is empty of data) and has 
little effect on the fit (not surprisingly, leaving out this correction increases \ 2 )- 

The line "ignore F — a p correlation" in Table 1 shows the change in the fitted parameters if 
we naively use the given noise estimates for weighting without accounting for the fact that there 
is a correlation between the flux estimate and the noise amplitude estimate for each pixel. Figure 
31 shows that the bias is a fairly constant 3-5% increase in the flux power. The difference is not 
actually caused by the change in weighting used in the power spectrum estimation - instead, the 
power is biased high because the estimation of the mean that each chunk of spectrum is divided 
by is biased low because low-flux pixels have smaller noise estimates. Ignoring this effect does not 
change our fit results. Normally we base our estimation of the amount of the noise that is due to 
quasar flux on the separate estimates we have from the spectral reduction pipeline for the flux, sky, 
and read- noise contributions, however these estimates do not add up to the total noise reported 
by the pipeline. If we rescale the individual numbers to make them consistent with the total (not 
necessarily the correct thing to do) we see that the fit results are not changed significantly (the 
line "rescale <r c " - we use a c to refer to noise computed using the separate flux, sky, and read-noise 
estimates), although the power does change by as much as 3% (Fig. 32 - this difference would be 
a bit larger if we did not directly measure and correct for the cross-correlation between the noise 
amplitude and flux). 

Finally, our power spectrum extraction code has a bias (partially related to division by the 
mean of each chunk of spectrum), that we correct for by dividing the result by Equation 11. 
Removing this correction decreases the estimated n c g by about 1/2 a and increases the amplitude 
by a similar amount. The combined change is actually quite significant because it is transverse to 
the degeneracy direction for these parameters. 
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Fig. 28. — Ratio of -Pio4i,ii85 computed without dividing by the continuum estimate to the standard 
case. 
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Fig. 29. — Ratio of -Pio4i,ii85 computed without the large constant added to the weight matrices 
to make the results less sensitive to the mean of the chunks, relative to the standard case. 
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Fig. 30. — Ratio of -Pio4i,ii85 computed with no correction for the offset between the defined center 
of each redshift bin and the center of weight of the data to the standard case. 
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Fig. 31. — Ratio of -Pio4i,ii85 computed without accounting for the correlation between the noise 
amplitude and the observed flux to the standard case. 
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Fig. 32. — Ratio of P1041 1185 computed using an alternative estimate of the fraction of the noise 
that is due to photon counting noise associated with flux from the quasar (see text) to the standard 
case. 



-65- 



To summarize, most of the effects described above are small relative to the statistical errors 
on the final estimated parameters. We understand the cases where the difference is significant, 
and expect that our standard method will be much more accurate than the difference between it 
and the variants (we show these variants to help the reader better understand our measurement). 
These tests give us confidence that the final results are very robust to small changes in the analysis 
pipeline. 

How sensitive are these conclusions to our assumptions about the nuisance parameters, F, 
T1.4, and 7, i.e., if these constraints improve in the future, will we need to worry more about 
systematic errors in Pp(k,z)'! We investigate this by first fixing all the parameters in the fit 
(including removing the noise amplitude uncertainty, resolution uncertainty, and freedom in the 
Silll correction), so the only uncertainty is on Pp(k,z). Table f (the "fixed nuisance parameters" 
line) shows that the error on the amplitude improves dramatically, by a factor of 5. The error on 
n e e improves by a factor of 2. So in principle the amplitude error can be improved a lot relative 
to potential systematic errors, and n e ff improved as well (see also Mandelbaum et al. 2003). The 
next line ("optimistic F"), where we assume the HIRES constraint on F is improved by a factor 
of 5, and the SDSS constraint by a factor of 2, shows that A 2 is substantially degenerate with F 
(as expected), but n e g is less degenerate. Improving the constraints on T1.4 and 7 by factors of 
5, in addition to the improved constraints on F, leads to little further improvement. Finally, for 
comparison, we tried simply reducing the errors on Pp(k,z) by a factor of \/3, and found that 
the error on n e fr decreases by almost the same factor (1.6), but the error on A 2 decreases less (a 
factor of 1.2). SDSS will collect a factor of ~ 3 more data than we have in the present sample. 
We conclude that the error on n e g- can easily be reduced by simply gathering more data, while 
improvements on A 2 can be made by improving the errors on F. 
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Table 1 — Continued 

Variant a 5 In A 2 cr lnA 2 5n cS cr„ cff Ax 2 b X 2 c 

P F (k,z) errors divided by ^3 -0.051 0.081 -0.000 0.015 



Note. — z p = 2.6, k p = 0.009 (kms" 1 )" 1 . 

a The meaning of each variant is explained in §4.3. 

b A^ 2 of the fitted parameters relative to the standard parameters, using the errors from the 
variant fit. 

C X 2 for the fit (essentially unrelated to Ax 2 ). 
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4.4. Subsamples of the Data 

Another way to test for systematic errors is to search for internal discrepancies between the 
different subsamples of the same data. Of course, there are only a finite number of possible subsam- 
ples we can try, so this test cannot be fully exhaustive. In addition, with many such tests performed 
one must worry that some will give an apparently statistically significant deviation just by random 
chance. Table 2 shows results of splitting the data into roughly equal weight subsamples, defined 
by various properties of the spectra that, at least at first glance, should not be correlated with the 
measured power. In practice, we rank the spectra by the property of interest and split the sample 
into halves by requiring that the Gaussian errors on the k = 0.007 (kms -1 ) -1 point are equal 
for the two halves (the bootstrap errors will not be precisely equal). We list the probability of 
obtaining x 2 greater than the value computed by differencing the power spectra (these differences 
include the different background subtraction computed using eqn. 13 for different noise levels). We 
also list the fitting parameter results for each subsample, and give the probability for obtaining the 
observed level of difference between the fits. Because these subsamples are basically independent, 
deviations within the error bars are expected and are not an indication of systematic errors. We 
describe these subsample splits below. 

The power we measure should be independent of the rest frame region of the quasar continuum 
in which it is measured. Figure 33 shows P1041, 1113 (&, z) and -Pm3,ii85(&, z) to test this expectation. 
The results look pretty similar, but to compare them quantitatively, we compute x 2 = (P< — 
P>) T (C< + C>) _1 (P< — P>), finding x 2 = 111-0 for 108 points. The agreement appears perfect. 
To compare the two in a different way, we perform separate fits of the linear mass power spectrum 
parameters A 2 and n e g to Pp(k, z) computed from Pio4i,iii3(&, z) and -P1113, 1185(^5 z). The results, 
given in the first line of Table 2, are consistent within the expected errors. This test provides some 
evidence that power from continuum fluctuations is not an important contribution to the total, 
beyond what we would expect from looking at the red side of the Lya emission line. It is possible 
that the two halves of the forest could have significant extra continuum power, but if they do it 
has to be the same in each half. 

We compute the weighted mean of the rms noise for each chunk of spectrum as we use it to 
estimate the power spectrum. A split based simply on this noise level, illustrated in Figure 34, 
produces a small but unambiguously significant discrepancy in the raw measurement of -Pio4i,ii85> 
X 2 = 185, though the fit parameters agree within their errors (Table 2, line 3). This discrepancy in 
power is the motivation for, and is largely removed by, our noise-dependent background subtraction 
procedure defined by equation (13). Figure 35 shows the power -Pi268,i380 that is used for background 
subtraction, again subsampled by noise level. There is a clear difference in power, which is not 
isolated to a few wavenumbers or redshifts. Once the -Pi268,i380 power is subtracted according to 
equation (13), we obtain Pp estimates and corresponding fit parameters from the high and low 
noise subsamples that agree within the errors (Table 2, line 2). Since the fit parameters agree even 
without noise-dependent background subtraction, it appears that the discrepancy in raw Pio4i,ii85 
power does not mimic a change in cosmological parameters, and our ultimate conclusions would 
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Table 2. Comparison between subsamples of the data 
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Probabilities may not be fully reliable because we have not demonstrated that x 2 is properly 
distributed. 

b The subsample fit results cannot be combined to produce the result of the fit to the full data 
set because the underlying nuisance parameters were not required to be the same. 

c The "noise (raw)" line shows the comparison without accounting for the noise dependence of 
the background. 



0.001 



0.01 



k [(krxi/s)- 1 ] 



Fig. 33. — Comparison of Pio4i,iii3(&, z) (connected by dotted lines) to i 2 ui3,ii85(&, z) (connected 
by solid lines, and shifted slightly to the right). The different redshifts have been shifted vertically 
by arbitrary amounts on this logarithmic plot (z increases from bottom to top). 
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Fig. 34. — -Pio4i,ii85 split by noise in the Lya forest region. The dotted line connects the low noise 
results, while the high noise results are offset slightly to the right. The results at different redshifts 
have been shifted vertically by arbitrary amounts (z increases from bottom to top). 
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Fig. 35. — -Pi268,i380 split by noise in the same region. The dotted line connects the low noise 
results, while the high noise results are offset slightly to the right. The results at different redshifts 
have been shifted vertically by arbitrary amounts. 
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therefore not change even if we did not implement it. Nonetheless, the origin of the difference 
remains somewhat mysterious, since we went to great effort to estimate noise correctly. 

Two more splits that yield discrepant Pio4i,ii85 but show no sign of trouble after the noise- 
dependent background subtraction are based on the ratio of the mean sky flux to the mean quasar 
flux and on a w — a c , the difference between the pipeline estimate of the noise and the sum of our 
estimates of the quasar flux, sky flux, and read-noise components of the noise. We are not sure 
what a w — a c means, since we do not understand the source of noise misestimation in the standard 
pipeline. Even without noise dependent background subtraction, the fit results did not differ 
significantly in either of these cases. They are almost surely symptoms of the same noise-related 
problem discussed above. 

The split based on read-noise in the spectra shows good agreement between the Pp{k,z) 
measurements, even without noise dependent background subtraction as does a split based on 
how well the mean continuum matches the quasar spectrum outside the Lya forest, quantified by 
computing x 2 /^ f° r the difference between the continuum and spectrum ("cont. x 2 /^" i n Table 
2). Several other splits that show little or no sign of trouble are based on: the mean value of x 2 ' l v 
computed for each pixel when combining exposures (this was the comparison that motivated our 
spectrum-by-spectrum noise re-estimation), the mean resolution, the movement of the spectrum 
relative to the detector pixel grid during the observation ( "flexure" ) , the alignment of the pixels in 
the different exposures for the same spectrum (closely related to flexure), the error on the overall 
normalization of the spectrum, A q (see eqn. 2, this error is set by a combination of the noise level 
outside the forest and the length of spectrum observed outside the forest), and the error on the 
means computed for the forest chunks (differences at fixed z are related to the length of the chunk 
and the noise in the forest). 

Overall, the agreement between our subsamples is excellent, both for the Pp(k,z) results and 
the fit results. In some cases this agreement relies on the noise-dependent background subtraction, 
which we would like to understand better (in no case does the fit agreement rely on this). 

4.5. Continuum Power 

The power in the mean continuum, for the 4 different rest frame regions identified in Figure 
2, is shown in Figures 36a-d, relative to the Lya forest power (the mean continuum power was 
measured by replacing the quasar flux in each pixel by the mean continuum level at that pixel). 
The mean continuum is very well behaved over the k range that we use (0.0013 — 0.02 (kms -1 ) -1 ), 
but its fluctuations quickly become significant at k < 0.001 (kms -1 ) -1 . What little power the 
mean continuum shows in our chosen k range should be removed when we divide the spectra by 
the continuum; it is only fluctuations around the mean that matter. 

We summarize our strong, but maybe not airtight, argument for believing that continuum 
fluctuations are not corrupting our measurement as follows: 
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0.001 0.01 



k [(krn/s)- 1 ] 

Fig. 36. — Power in the mean continuum relative to the Lya forest power, for various rest wave- 
length intervals: (a) 1041 < A rcst < 1185A, (b) 1268 < A rcst < 1380A, (c) 1409 < A rcst < 1523A, (d) 
1558 < A rcst < 1774A. The error bars show the fractional error on -Pio4i,ii85 (without diagonalizing 
the window matrix because the diagonalization works poorly at the lowest A:s that we show). 
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• The power in the mean continuum is small. 

• The results for the Pp{k,z) measurement in two halves of the forest region, Pio4i,ni3 and 
^1113,1185, agree. 

• The power in the background regions, -Pi268,i380 and -Pi409,i523, agree at the level we care 
about in the low-noise data, as long as BALs (which mostly affect the latter region) are 
removed. 

• Division by preliminary 13 eigenvector PC A estimates of the continua (i.e., including fluctu- 
ations around the mean) does not change the results. 

To be quantitatively important despite these arguments, the power in quasar-to-quasar continuum 
fluctuations in the forest must be substantially larger than the power in the mean continuum 
itself, the continuum fluctuations in the forest must be substantially different from those in the 
background regions (despite those regions being similar to each other and the two halves of the 
forest being similar to each other), and our PC A analysis must be substantially flawed. Further 
study is warranted, but a big effect seems unlikely. 

4.6. Comparison with Past Measurements 

There are three Pp(k, z) measurements already in the literature, McDonald et al. (2000), Croft 
et al. (2002), and Kim et al. (2004), all using at least some high resolution data. Each uses its 
own set of redshift bins, so to compare we need a way to interpolate our results to these redshifts. 
We do this by performing our standard cosmological fit to all of the data (at first - later we will 
remove some of the past results). This gives us a set of best fit model parameters that can be used 
to compute the power at any k and z. Within the range of k where we have SDSS measurements, 
the fit is always dominated by the SDSS points. The fitted curves always match the SDSS results 
to much better than the size of the errors on the past results, meaning that, for the the purpose 
of comparison to the past results, the curves are simply a faithful interpolation between the SDSS 
points. At k > 0.02(kms -1 ) -1 , the fit is effectively a weighted average of the past results, although 
the constraint that it must match SDSS at lower k has some influence (our simulation predictions 
do not allow for sharp features in Pp(k,z).) 

We first perform a fit to all the data with k < 0.05(kms -1 ) -1 and z > 2.1, finding an 
atrociously bad \ 2 = 392 for ~ 238 dof. Removing McDonald et al. (2000), reduces \ 2 by 53.4 
(for 39 data points), removing Croft et al. (2002) reduces x 2 by 85.2 (for 65 points), and removing 
Kim et al. (2004) reduces x 2 by 123.3 (28 points). Clearly there is gross disagreement between Kim 
et al. (2004) and the other results. Figure 37 shows the Kim et al. (2004) points at z = 2.18 and 
z = 2.58 (from their Table 5) along with the fit prediction for them. Note that we include Silll 
contamination in the model as described in §4.2, so the model curves are not perfectly smooth. 
We see large discrepancies, as we expect from the bad \ 2 ■ The point at k = 0.0012( kms -1 ) -1 , 
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Fig. 37. — The black, solid line and open squares (blue, dotted line and crosses) show the fit 
prediction and measured points from Kim et al. (2004) at z = 2.18 (2.58). Red (green) error bars 
show our SDSS measurement at z = 2.2 (2.6). 
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z = 2.58 is 5.9a below the prediction (as well as any reasonable extrapolation of the other Kim 
et al. (2004) points), and the points at increasingly high k are generally too high (at the highest k, 
this is a reflection of disagreement with the other high resolution data, but actually the agreement 
is not much better if we only include SDSS in the fit, because no model can fit the highest k SDSS 
points and then climb to match the higher k Kim et al. (2004) points). To reassure the reader that 
we are not playing games with the fitted curves, we also plot the SDSS points at z = 2.2 and 2.6. 

Since the Kim et al. (2004) results clearly have some problem, unless the other three measure- 
ments are all wrong (we will see that, with one exception, the other three agree with each other), 
we eliminate them from the rest of the comparison. A fit to SDSS, McDonald et al. (2000), and 
Croft et al. (2002) gives \ 2 = 269 for ~ 210 dof (still a bad fit). Removing McDonald et al. (2000), 
reduces x 2 by 44.3 (39 points), while removing Croft et al. (2002) reduces x 2 by 97.2 (for 65 points 
this reduction would occur by chance only 0.6% of the time). Figure 38 shows the Croft et al. 
(2002) points, along with the fit prediction for them. The agreement is actually very good for 4 of 
the 5 redshift bins, while the z = 2.47 points are obviously out of place (these 13 points increase 
X 2 by 54). Figure 39 shows the McDonald et al. (2000) points, along with the fit prediction for 
them (for this fit we removed the z = 2.47 Croft et al. (2002) points). The agreement is good, 
with the agreement at z = 2.41 disfavoring the anomalous Croft et al. (2002) z = 2.47 points; 
further investigation by R. Croft (private communication) does not reveal any obvious error in this 
redshift bin that would explain the anomaly. Note that the agreement of McDonald et al. (2000) 
and Croft et al. (2002) at high k adds weight to the idea that something is seriously wrong with 
Kim et al. (2004). Kim et al. (2004) show some comparisons with past results, and claim they 
agree, but these comparisons used custom redshift bins (i.e., not the bins in their table), and were 
not high precision (for example they compare the Croft et al. (2002) points at z = 2.13 to a bin 
with z = 2.04, so evolution cancels some of the amplitude offset, and they call the apparent ~ 50% 
difference at k ~ 0.04( kms -1 ) -1 a "slight" excess). 



5. Final Results Table and Directions for Use 

Table 3 gives the primary power spectrum results. The columns are: z, the redshift of the bin; 
k, the wavenumber of the bin; Pp(k, z), our final Lya forest power spectrum result (along with the 
square roots of the diagonal elements of the error covariance matrix); -P no i S e) the noise power that 
was subtracted from each bin; and -Pbackground , the background that was subtracted from each bin 
(-Pi268,i380 adjusted according to the amount of noise in the forest, eqs. 3, 4, and 13). -Pnoise is just 
the noise subtracted from -Pio4i,io85 (a roughly comparable amount of noise was subtracted from 
the background, so to some degree these cancel in the final result). The table and the covariance 
matrix of the errors are available at http://feynman.princeton.edu/~pmcdonal/LyaF/sdss.html. 
The covariance matrix must be used in any serious quantitative fitting. When using this table to 
constrain models, the following allowances should be made for residual systematic uncertainties: 



-81- 




Fig. 38. — Measured points and fit prediction for the Croft et al. (2002) results. From bottom 
to top (roughly) — z=2.13: black, solid line, open square; z=2.47: blue, dotted line, 4-point star 
(cross); z=2.74: cyan, dashed line, filled square; z=3.03: green, long-dashed line, open triangle; 
z=3.51: magenta, dot-dashed line, 3-point star. 
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Fig. 39. — Measured points and fit prediction for the McDonald et al. (2000) results. ^From bottom 
to top (roughly) — z=2.41: black, solid line, open square; z=3.00: blue, dotted line, 4-point star 
(cross); z=3.89: cyan, dashed line, filled square; 
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• Allow ±5% rms error on the noise-power amplitude at each redshift. We do not have any 
reason to think the error is really this large, but, considering the complications related to the 
noise, we think it is prudent to include it. Operationally, we suggest subtracting fiP U oise(k, z^) 
from Zi), where /j are free parameters in the fit (one for each redshift bin), and adding 
£^/0.05) 2 to x 2 . 

• Allow ±(7kms~ 1 ) 2 rms overall error on the resolution variance (i.e., the square of the rms 
width of the Gaussian resolution kernel) . This is the expected size of the uncertainty due to 
flexure in the detector, although Figure 5 suggests that it may actually be smaller. Specifically, 
multiply Pp(k,z) by exp(afc 2 ), with a a free parameter in the fit, and add [a/(7kms -1 ) 2 ] 2 
to x 2 - 

• Silll-Lya cross-correlation must be accounted for. We have suggested a simple method - 
assume the cross-correlation has the same form as the Lya-Lya auto-correlation up to an 
amplitude that is a free parameter, and possibly include freedom in the correlation width 
and/or redshift evolution of the amplitude - but others could be devised (e.g., including Silll 
in the simulated spectra through a parameterized semi-analytic model) . 
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Table3. P F (k, z) Results 



Pj?{k^ z} noise -^background 



2.2 


0.00141 


18.09 ± 1.74 


6.20 


3.29 


2.2 


0.00178 


17.55 ± 1.34 


6.07 


2.71 


2.2 


0.00224 


19.05 ± 1.21 


6.20 


2.46 


2.2 


0.00282 


18.93 ± 1.03 


6.23 


2.45 


2.2 


0.00355 


15.80 ± 0.74 


5.92 


1.97 


2.2 


0.00447 


12.68 ±0.60 


5.57 


1.18 


2.2 


0.00562 


14.04 ±0.53 


5.85 


1.04 


2.2 


0.00708 


11.09 ± 0.45 


5.87 


1.17 


2.2 


0.00891 


9.38 ± 0.39 


6.06 


1.20 


2.2 


0.01122 


8.09 ± 0.38 


6.87 


1.40 


2.2 


0.01413 


6.99 ± 0.34 


8.48 


1.32 


2.2 


0.01778 


4.69 ± 0.35 


10.52 


0.87 


2.4 


0.00141 


21.52 ± 1.91 


5.70 


3.72 


2.4 


0.00178 


23.66 ± 2.09 


5.65 


2.63 


2.4 


0.00224 


23.57 ± 1.39 


5.65 


2.45 


2.4 


0.00282 


22.25 ± 1.24 


5.58 


2.47 


2.4 


0.00355 


18.65 ±0.89 


5.31 


1.82 


2.4 


0.00447 


15.74 ±0.65 


5.15 


1.18 


2.4 


0.00562 


18.07 ±0.68 


5.43 


0.79 


2.4 


0.00708 


13.16 ±0.51 


5.32 


1.07 


2.4 


0.00891 


12.58 ±0.41 


5.71 


1.14 


2.4 


0.01122 


10.42 ±0.41 


6.27 


1.29 


2.4 


0.01413 


8.17 ±0.36 


7.32 


1.22 


2.4 


0.01778 


6.08 ±0.33 


9.37 


0.91 


2.6 


0.00141 


28.29 ±2.55 


6.78 


4.02 


2.6 


0.00178 


29.04 ± 1.85 


6.68 


3.18 


2.6 


0.00224 


32.13 ± 1.76 


6.76 


2.65 


2.6 


0.00282 


27.44 ± 1.39 


6.63 


2.41 


2.6 


0.00355 


25.06 ± 1.09 


6.52 


1.79 


2.6 


0.00447 


20.67 ±0.85 


6.40 


1.27 


2.6 


0.00562 


22.49 ± 0.72 


6.69 


1.05 


2.6 


0.00708 


17.19 ±0.60 


6.71 


1.17 


2.6 


0.00891 


15.40 ±0.51 


7.12 


1.05 



Table 3 — Continued 



-fnoise -^background 



O f* 

2.6 


0.01122 


1 o n r* i n A o 

13.25 ± 0.48 


*~7 n 1 

7.91 


1 O A 

1.24 


2.6 


n t\\ a i o 

0.01413 


i n or i n a i 

10.25 ± 0.41 


9.15 


1 OO 

1.22 


2.6 


0.01778 


8.43 ± 0.37 


1 1 (' (' 

11.66 


n n r 

0.95 


2.8 


0.00141 


37.25 ± 2.75 


f* O O 

6.83 


o nn 

2.99 


2.8 


n nm 'vo 

0.00178 


o ^7 r o i o on 

37.52 ± 2.20 


6.75 


o i n 

2.10 


2.8 


/-\ onn O A 

0.00224 


q o T A 1 1 on 

38.74 ± 1.80 


6.79 


1.84 


2.8 


0.00282 


o T 1 o i -i /i o 

37.12 ± 1.48 


6.78 


2.29 


o o 

2.8 


0.00355 


on 11 ii oo 

30.11 ± 1.23 


n c n 

6.60 


i r o 

1.52 


2.8 


n nn A A ^ 

0.00447 


or - ci i n no 

25.67 ± 0.92 


6.52 


1 OO 

1.22 


2.8 


0.00562 


or f7j i n o o 

25.74 ± 0.83 


6.73 


1.19 


2.8 


0.00708 


OO V A in C T 

22.54 ± 0.67 


6.95 


0.98 


2.8 


r\ on O n 1 

0.00891 


on 1 o i n (■ • o 

20.12 ± 0.62 


7.41 


1.11 


O O 

2.8 


n n i 1 oo 

0.01122 


i r on i n /i o 

15.89 ± 0.48 


8.06 


1 1 O 

1.18 


o a 

2.8 


0.01413 


i o n /i i n A o 

13.04 ± 0.42 


n o ^7 

9.37 


1 1 o 

1.13 


o o 

2.8 


0.01778 


c •> / ' • > i no/ 1 

9.63 ± 0.36 


11.64 


0.90 


3.0 


r\ f\r\-\ A 1 

0.00141 


/i n o <^ i o ^70 

46.36 ± 3.72 


7.76 


3.51 


3.0 


0.00178 


42.53 ± 2.87 


7.63 


2.74 


3.0 


0.00224 


47.66 ± 2.69 


7.73 


2.20 


3.0 


0.00282 


/i o on i o i n 

42.20 ± 2.19 


7.66 


2.84 


o n 

3.0 


r\ nn o r r 

0.00355 


o ri rin i i 'to 

36.99 ± 1.72 


7.51 


1 O 1 

1.81 


3.0 


0.00447 


29.47 ± 1.20 


7.34 


1.30 


3.0 


0.00562 


30.12 ± 1.07 


7.56 


1.33 


3.0 


0.00708 


24.30 ± 0.81 


7.63 


0.99 


3.0 


0.00891 


22.51 ±0.75 


8.15 


1.30 


3.0 


0.0H22 


18.75 ±0.66 


8.83 


1.30 


3.0 


0.01413 


14.33 ±0.52 


9.89 


1.38 


3.0 


0.01778 


11.26 ±0.47 


11.90 


0.91 


3.2 


0.00141 


54.73 ± 4.97 


9.57 


5.07 


3.2 


0.00178 


49.72 ±4.12 


9.44 


3.73 


3.2 


0.00224 


52.86 ± 3.29 


9.46 


3.01 


3.2 


0.00282 


48.44 ± 2.58 


9.38 


2.59 


3.2 


0.00355 


44.01 ± 2.33 


9.28 


2.71 


3.2 


0.00447 


35.12 ± 1.54 


8.95 


1.36 
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Table 3 — Continued 



-fnoise -^background 



o o 

3.2 


n nn f / ' o 

0.00562 


O A r T 1 1 A 1 

34.57 ± 1.41 


9.15 


1 o o 

1.28 


3.2 


0.00708 


1 1 A 1 1 1 A 

31.14 ± 1.19 


n a n 

9.40 


1 on 

1.29 


o o 

3.2 


n AAOA1 

0.00891 


26.9b ± 0.95 


i\ on 

9.80 


1 A O 

1.48 


o o 

6.2 


A A1 1 OO 

0.01122 


OA 01 I A OO 

22.21 ± 0.83 


1 n /in 

10.49 


1 on 

1.80 


o o 

6.2 


n ni a i o 

0.01413 


1 O O T 1 n Tn 

18.37 ± 0.70 


11 T A 

11.74 


1 o o 

1.38 


3.2 


0.01778 


15. 12 ± 0.66 


14.07 


1.18 


3.4 


0.00141 


p* ir» i r* o p* 

56.42 ± 5.85 


11.12 


4.08 


O A 

3.4 


r\ r\r\ -i TO 

0.00178 


TP" Tr i r o o 

75.75 ± 5.33 


1 1 OO 

11.32 


O A 1 

2.41 


O A 

3.4 


n aaoa A 

0.00224 


r*/ 1 Tn i o ot 

56.79 ± 3.87 


10.97 


2.09 


O A 

3.4 


n AAOOO 

0.00282 


f o 1 n i o /in 

58.40 ± 3.43 


11.04 


3.16 


O A 

3.4 


n nn or r 

0.00355 


52.56 ± 2.85 


10.96 


2.62 


O A 

3.4 


n r\r\ a A*7 

0.00447 


43.43 ± 2.21 


10.76 


2.00 


O A 

3.4 


r\ r\r\ — / < o 

0.00562 


A~\ 1 1 TO 

41.67 ± 1.73 


1 n nn 

10.99 


o n A 

2.04 


3.4 


n nn ^7A o 

0.00708 


O T O/ 1 1 1 /I O 

37.36 ± 1.43 


1 1 oo 

11.28 


1 f P* 

1.55 


O A 

3.4 


0.00891 


o o r 1 T i i in 

32.57 ± 1.19 


11.87 


1 TT 

1.77 


3.4 


n a i i ■ \c\ 

0.01122 


O O P* 1 1 1 "1 P* 

28.51 ± 1.15 


13.06 


2.00 


O A 

3.4 


0.01413 


oo i n oo 

22.28 ± 0.88 


14.63 


1.63 


3.4 


0.01778 


18.01 ± 0.79 


17.80 


1.23 


3.6 


0.00141 


to ,1 r* i o o o 

79.46 ± 8.33 


15.11 


2.25 


3.6 


r\ nn i to 

0.00178 


Of in I o OO 

85.12 ± 8.28 


i A nn 

14.90 


o OT 

2.87 


3.6 


0.00224 


75.03 ± 5.88 


14.87 


3.28 


3.6 


0.00282 


66.15 ±4.98 


14.51 


3.30 


3.6 


0.00355 


66.32 ±4.08 


14.59 


2.26 


3.6 


0.00447 


55.66 ± 3.36 


14.22 


1.33 


3.6 


0.00562 


49.51 ± 2.72 


14.17 


1.28 


3.6 


0.00708 


43.77 ± 2.15 


14.41 


1.62 


3.6 


0.00891 


40.20 ± 1.93 


15.09 


1.98 


3.6 


0.01122 


32.04 ± 1.63 


15.72 


1.61 


3.6 


0.01413 


25.82 ± 1.31 


17.26 


1.25 


3.6 


0.01778 


21.49 ± 1.23 


21.11 


1.51 


3.8 


0.00141 


118.61 ± 15.47 


22.58 


6.89 


3.8 


0.00178 


61.52 ±9.80 


20.93 


6.40 


3.8 


0.00224 


77.29 ± 7.09 


20.91 


5.07 
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Table 3 — Continued 



z 


k 


P F (k,z) 


1 noise 


-^background 


3.8 


0.00282 


71.78 ± 7.42 


20.76 


2.20 


3.8 


0.00355 


77.49 ± 5.65 


21.00 


2.51 


3.8 


0.00447 


59.12 ±4.21 


20.37 


2.37 


3.8 


0.00562 


57.53 ± 3.72 


20.78 


2.61 


3.8 


0.00708 


56.25 ± 3.45 


21.63 


2.79 


3.8 


0.00891 


42.46 ± 2.34 


21.69 


2.57 


3.8 


0.01122 


36.93 ± 2.25 


23.43 


2.69 


3.8 


0.01413 


29.52 ±2.25 


26.14 


2.47 


3.8 


0.01778 


27.72 ± 1.85 


33.51 


2.15 



Note. — k has units (kms ) 1 , power spectra have units kms . The error covariance matrix 
must be used for any quantitative fitting. 
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6. Conclusions 

We have analyzed a sample of 3035 quasar spectra measured by SDSS and covering the redshift 
range 2 < z < 4. This data set is almost two orders of magnitude larger than previously available 
data sets. We have focused on the flux power spectrum in the redshift range 2.1 < z < 3.9 and for 
modes 0.0013 (kms -1 ) -1 < k < 0.02 (kms -1 ) -1 . The extraordinary size of the data sample leads 
to an order of magnitude reduction in errors compared to previous analyses. Consequently, to do 
justice to this data set requires a much more careful analysis than was needed in the past. To this 
end we have developed a new analysis pipeline using quadratic power spectrum estimation with 
near optimal performance. We applied this analysis to realistic mock spectra and demonstrated 
(after several tweaks) that the method performs as expected. We emphasize that realistic mock 
spectra are essential if one is to trust the results at the level of precision allowed by this data 
set. Our error estimation is based on bootstrap resampling, which works well here because the 
individual quasars are independent of each other. The errors were tested against mock spectra and 
found to be accurate. We also compared the bootstrap errors to Gaussian errors, finding them to 
be in general less than 20% higher than Gaussian. 

Given the small errors on the recovered flux power spectrum the required control of systematic 
effects must be improved correspondingly as well. A significant part of this paper is devoted to this 
issue. We find several sources of contamination present in the data and develop methods to remove 
them. Metal absorption for metals with rest wavelength transitions above ~ 1300A, uncertainties 
in sky subtraction, and calibration errors can be subtracted essentially exactly by measuring the 
power on the red side of Lya forest using the same observed wavelength range. We search for 
a contribution from metals with transitions close to Ly-a using a correlation function analysis, 
assuming that they are correlated with hydrogen. We find clear evidence of Silll contamination 
and develop a simple and effective scheme to remove it. This procedure improves the \ 2 of the fit 
from 194 to 129 for ~ 104 degrees of freedom and is thus necessary for a satisfactory fit. We find no 
evidence of any other metal line contribution to the background subtracted flux power spectrum. 

We reduce any contribution of the continuum to the flux power spectrum by dividing each 
spectrum by the mean quasar continuum. If contributions from quasar-to-quasar continuum dif- 
ferences are similar in different regions of the spectrum, then our subtraction of power from the 
red side of Lya, as described above, should remove them. Several tests suggest that any residual 
contributions from continuum fluctuations are negligible. First, we measure the power in the mean 
continuum in several rest frame regions, finding it to be always small relative to our error bars, so 
power in quasar-to-quasar fluctuations has to be larger than power in the mean continuum itself to 
be significant. Second, measurements of the background in several rest frame regions place upper 
limits on the fluctuations in those regions. Third, a split of the Lya forest region into two halves 
reveals no evidence that residual continuum fluctuations differ from one half to the other. Finally, 
estimating the continuum quasar-by-quasar using a principal component analysis does not change 
the power spectrum results significantly. 
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In section §4 we perform a series of tests to verify the robustness of the analysis against 
several modifications of the standard procedure and splits of the data. This reveals an interesting 
correlation between the power in the red side and the average noise (and some other properties of 
the spectra that correlate with noise, like the amplitude of the sky flux relative to the quasar flux). 
While we do not have a detailed explanation for this effect, we are able to remove it by modifying 
the standard procedure to include this correlation. From our full battery of tests, we conclude 
that systematic effects in the power spectrum measurement are not likely to significantly affect the 
results of cosmological fitting (i.e., it is likely that some effects remain formally significant relative 
to the errors on Pp(fc, z), but the shape of these systematic errors does not seem to correspond to 
a change in the cosmological parameters). This conclusion is further confirmed by the analysis of 
different subsets, which do not reveal any systematic deviations from those expected statistically. 

In this paper we limit ourselves to the analysis of the flux power spectrum, without attempting 
to compare it to cosmological models. The results of this paper should thus be fairly noncontro- 
versial and can be used by others who wish to perform their own cosmological analysis. Our own 
analysis will be presented in a separate publication, as will the cosmological implications that follow 
from it. We note that the expected error on the linear rms amplitude of fluctuations is ~ 5% and 
on the slope is ~ 0.024, both at the pivot point k = 0.009 (kms -1 ) -1 . This should be compared to 
10% error on the amplitude and 0.04 on the slope from the WMAP data at k = 0.05Mpc _1 (Spergel 
et al. 2003). This data set provides very tight constraints on the amplitude and slope of the matter 
power spectrum. Many additional analyses can be performed using this data set, among them 
the mean flux evolution, cross-correlations between close pairs, and a bispectrum measurement. 
These will provide a wealth of additional information both on cosmology and on the state of the 
intergalactic medium at 2 < z < 4, and they will allow us to test the basic picture of the Lya forest 
that has emerged over the last decade. We believe that the unprecedented size of this data set will 
revolutionize our understanding of the high redshift universe; this work is merely a first step in this 
endeavor. 
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